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ABSTRACT 

We present a comprehensive description of the population synthesis code StarTrack. The original 
code has been significantly modified and updated. Special emphasis is placed here on processes leading 
to the formation and further evolution of compact objects (white dwarfs, neutron stars, and black holes). 
Both single and binary star populations are considered. The code now incorporates detailed calculations 
of all mass-transfer phases, a full implementation of orbital evolution due to tides, as well as the most 
recent estimates of magnetic braking. This updated version of StarTrack can be used for a wide 
variety of problems, with relevance to many current and planned observatories, e.g., studies of X-ray 
binaries (Chandra, XMM-Newton), gravitational radiation sources (LIGO, LISA), and gamma-ray burst 
progenitors (HETE-II, Swift). The code has already been used in studies of Galactic and extra-galactic 
X-ray binary populations, black holes in young star clusters, Type la supernova progenitors, and double 
compact object populations. Here we describe in detail the input physics, we present the code calibration 
and tests, and we outline our current studies in the context of X-ray binary populations. 

Subject headings: binaries: close — stars: evolution — stars: white dwarfs, neutron — black hole physics 
— X-rays: binaries 



1. INTRODUCTION 

The StarTrack population synthesis code was initially 
developed for the study of double compact object merg- 
ers in the context of gamma-ray burst (GRB) progeni- 
tors (Bclczynski, Bulik & Rudak 2002b) and gravitational 
radiation inspiral sources (Belczynski, Kalogera & Bulik 
2002c, hereafter BKB02). StarTrack has undergone ma- 
jor updates and revisions in the last few years. With this 
code we are able to evolve isolated (not dynamically in- 
teracting) single stars and binaries for a wide range of ini- 
tial conditions. The input physics incorporates our latest 
knowledge of processes governing stellar evolution, while 
the most uncertain aspects are parameterized to allow for 
systematic error analysis. During the code development, 
special emphasis was placed on the compact object popula- 
tions: white dwarfs (WDs), neutron stars (NSs), and black 
holes (BHs). The input physics currently includes all ma- 
jor processes important for the formation and evolution of 
compact objects. Among other things we have developed 
fast procedures to treat and diagnose various types of mass 
transfer episodes (including phases of thermal timescale 
and dynamically unstable mass transfer leading to com- 
mon envelopes). We also compute tidal effects on orbital 
evolution, angular momentum losses due to magnetic brak- 
ing and gravitational radiation, as well as mass loss from 
stellar winds and during mass transfer phases. Rejuvena- 
tion of binary components is taken into account. The full 
orbital evolution of binaries is also computed, including 



angular momentum and mass loss. Supcrnovae (SNe) and 
compact object formation are also treated in detail. 

The new version of StarTrack presented here has al- 
ready been tested and used in many applications. Belczyn- 
ski & Taam (2004a) studied the formation of ultrashort pe- 
riod X-ray binaries and they also demonstrated that the 
faint X-ray Galactic Center population can neither be ex- 
plained by quiescent NS/BH transients nor by hard/faint 
wind-fed sources (Belczynski & Taam 2004b). Belczynski, 
Sadowski & Rasio (2004b) and Belczynski et al. (2006) de- 
veloped a comprehensive description of young BH popula- 
tions, which can also provide realistic initial conditions for 
the dynamical modeling of BHs in star clusters. Belczynski 
et al. (2004a) derived for the first time a synthetic X-ray 
luminosity function which agrees with Chandra observa- 
tions of NGC 1659, and Scpinsky, Kalogera, & Belczynski 
(2005) explored the numbers and spatial distribution of 
X-ray binaries formed in young star clusters. Belczynski, 
Bulik & Ruiter (2005b) tested different models of Type la 
SN progenitors, arriving at the conclusion that the double 
degenerate scenario most easily reproduces the observed 
delay times between star formation and Type la SNe. Bel- 
czynski et al. (2005a) used StarTrack to study the grav- 
itational radiation signal from the Galactic population of 
double WDs. Nutzman et al. (2004), O'Shaughnessy et al. 
(2005a, b,c), and Ihm, Kalogera, & Belczynski (2005) stud- 
ied binary compact object populations and derived merger 
rates and detection rates by ground-based interferometers; 
they also examined BH spin magnitudes and studied the 
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eccentricities of double neutron stars. StarTrack was also 
incorporated into a simple stellar dynamics code, allow- 
ing the study of the effects of dynamical interactions on 
binary populations in dense star clusters. In that form it 
has been used for the study of binary fractions in globu- 
lar clusters (Ivanova et al. 2005) and an investigation of 
intermediate-mass BHs in clusters and their connection to 
ultra- luminous X-ray sources (Blecha et al. 2005). 

Among other things StarTrack has been adapted for 
the study of accretion powered X-ray binaries (XRBs). In 
forthcoming papers we will present the synthetic popula- 
tions of XRBs formed in different stellar environments. We 
will start with young starburst galaxies, and move on to 
spiral and, eventually, old elliptical galaxies. In the next 
stage it will be possible to compare the models with rapidly 
improving observations of various X-ray point source pop- 
ulations. This will offer a new perspective to the study 
of several uncertain aspects of binary evolution leading to 
the formation of XRBs. It may also result in an indepen- 
dent diagnostic of star formation rates for nearby galax- 
ies, since both the numbers and properties of XRBs are 
directly connected to the star formation history (see e.g., 
Grimm, Gilfanov & Sunyaev 2003; Gilfanov 2004; Kim & 
Fabbiano 2004; Belczynski et al. 2006, in preparation). 

In this paper we provide a detailed description of the 
current version of StarTrack, and we present the results 
of a number of tests. We describe the implementations of 
single star evolution in §2, binary orbit evolution in §3, 
stellar wind mass loss/accretion in § 4, Roche lobe overflow 
calculations in § 5, spatial velocities in § 6, and the assumed 
distributions of initial parameters in § 7. In § 8 we discuss 
the validity of various input physics assumptions, and we 
compare StarTrack calculations with detailed evolution- 
ary models and with various observations. Section § 9 is 
dedicated to the discussion of X-ray binary modeling. In 
§ 10 wc conclude with a short summary. 

2. SINGLE STELLAR EVOLUTION 

In all subsequent sections we use units of Mq for mass, 
i?o for orbital separations and stellar radii, Myr for time, 
Lq for bolometric luminosity, unless specified otherwise. 
We use R and M to denote stellar radius and mass, while 
a, e represent the binary orbital parameters: semi-major 
axis and eccentricity, respectively. Index i = 1,2 is used 
to mark the binary components (or single stars for con- 
sistency), or to denote an accretor and a donor in mass 
transfer calculations: i = acc, don. Roche lobe parameters 
are indexed with "lob". The initially more massive (at 
Zero Age Main Sequence) binary component is referred to 
as primary, while its companion as secondary. 

2.1. Overview 

The evolution of single stars and non-interacting binary 
components have remained mostly unchanged since the 
last pubhshed description of the code (BKB02) and there- 
fore we only give a brief outline here. However, we do point 
out the new additions and reiterate the modifications to 
the original formulas which were used as the base for the 
implementation of single star evolution in StarTrack. 

To evolve single stars from the Zero Age Main Sequence 
(ZAMS) until remnant formation (WD, NS, BH, or a 
remnant-less supernova) we employ the analytic formulas 
of Hurley, Pols & Tout (2000). Each star is followed along 



an evolutionary track specific for its initial mass and metal- 
licity. Various wind mass loss rates that vary with the stel- 
lar evolutionary stage are incorporated into the code and 
their effect on stellar evolution is taken into account. Once 
the remnant is formed, we terminate the calculations but 
keep track of the numbers, properties and formation times 
of a given type of remnant. Additionally, for white dwarf 
remnants we take into account their subsequent luminosity 
evolution, and follow cooling tracks adopted from Hurley 
et al. (2000). 

2.2. Stellar types 

We follow Hurley et al. (2000) to denote different stages 
of stellar evolution with an integer Ki = l..n, where 

- Main Sequence (MS) M < 0.7 Mg 

1 - MS M > 0.7 Mq 

2 - Hertzsprung gap (HG) 

3 - Red Giant Branch (RG) 

4 - Core Helium Burning (CHeB) 

5 - Early Asymptotic Giant Branch (EAGB) 

6 - Thermally Pulsing AGB (TPAGB) 

7 - Helium Main Sequence (HeMS) 

8 - Helium Hertzsprung gap (HeHG) 

9 ~ Helium Giant Branch (HeGB) 

10 - Helium White Dwarf (He WD) 

11 - Carbon/Oxygen White Dwarf (CO WD) 

12 - Oxygen/Neon White Dwarf (ONe WD) 

13 - Neutron Star (NS) 

14 - Black Hole (BH) 

15 - massless remnant (after SN la explosion) 

16 - Hydrogen White Dwarf (H WD) 

17 - Hybrid White Dwarf (Hyb WD) 

In addition to the star types introduced and coded by 
the numbers Ki ~ 1...15 in the original Hurley et al. 
(2000) formulas, we have introduced two new stellar types 
Ki = 16, 17. Ki = 16 denotes a H-rich white dwarf. Only 
main sequence stars less massive than about 0.7 Mq can 
produce such a H-rich remnant through mass loss in a 
close binary system. These low-mass stars do not process 
a significant amount of hydrogen into helium in their cores 
(even in a Hubble time) and once their mass is stripped 
below the hydrogen burning hmit (close to ~ 0.08 M©) 
they become degenerate H-rich white dwarfs. These stars, 
although not frequently encountered in population synthe- 
sis, may become donors in the shortest-period interacting 
binaries. Ki = 17 denotes a hybrid white dwarf, with 
a carbon-oxygen-helium mixture in the core and a helium 
envelope. These objects are the remnants of naked Helium 
main sequence stars {Ki = 7) which are stripped of mass 
below 0.35 Mq during Roche lobe overflow (RLOF). At 
that point, thermonuclear reactions stop and the star be- 
comes degenerate (eg., Savonije, de Kool & van den Heuvel 
1986). 

2.3. Modifications 

Several major changes to the original Hurley et 
al. (2000) formulas have been implemented within 
StarTrack. 

2.3.1. Compact object masses 

The remnant masses of neutron stars and black holes are 
calculated in a different way than originally suggested by 
Hurley et al. (2000). In the present version of the code we 
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have further revised our prescription presented in Belczyn- 
ski et al. (2002c) to include the more recent calculations 
of FeNi core masses and allow for the possibility of NS for- 
mation through electron capture supernovae (ECS). White 
dwarfs masses are calculated with the original formulas of 
Hurley et al. (2000), although ONe WDs are formed in a 
slightly narrower range since we allow for ECS NS forma- 
tion (see below). 

We determine the mass of a NS/BH remnant using 
information on the final CO and FeNi core masses, com- 
bined with the knowledge of the pre-supernova mass of 
the star. For a given initial ZAMS mass, the final CO core 
mass is obtained from the original Hurley et al. (2000) 
formulas, while we use the models of Timmes, Woosley 
& Weaver (1996) to estimate final FeNi core mass. The 
results of Timmes et al. (1996, see their Fig. 2) show two 
distinctive FeNi core masses (we use models with the ad- 
dition of Si shell mass), below and above initial masses 
of Mzams ~ 18 — 19 M0. The dichotomy arises from dif- 
ferent carbon burning (convective versus radiative) the in 
pre-supernova stellar core. For higher mass progenitors 
(Afzams > 20 M0) there is a slow rise in a FeNi core mass 
that we approximate with a linear relation. The final FeNi 
core mass for a given CO core mass (Mco) is obtained 
from: 

-^'^FoNi = 

f 1.50 Mco < 4.82 (M,ams < 18) 

2.11 4.82 < Mco < 6.31(18 < M.^ms < 25) 

' 0.69Mco - 2.26 6.31 < Mco < 6.75(25 < A/^ams < 30) 

, 0.37Mco - 0.07 Mco > 6.75 (Af^ams > 30) 

where all masses are expressed in M0. The above depen- 
dence was obtained for solar metallicity models, but we use 
it for all compositions considered here {Z = 0.0001 — 0.3; 
the metallicity range covered by Hurley et al. 2000 in fits 
to the stellar evolution models). Timmes et al. (1996) 
results for metallicity are very similar for zero metallicity. 
For example, in the mass range most important for NS 
formation the core mass changes from ~ 1.50 M0 (solar 
metallicity) to ~ 1.58 M0 (zero metallicity). The differ- 
ences can be larger for BH progenitors, but then the exact 
core mass does not play such an important role on the fi- 
nal BH mass, as most BHs form through fall back or direct 
collapse (see below). 

The effects of material fallback (ejected initially in the 
SN explosion) during the star's final collapse are included. 
For the most massive stars we also allow for the possibility 
of a silent collapse (no supernova explosion) and direct BH 
formation. For solar metallicity and standard wind mass 
loss the compact object masses are obtained from 

r A/peNi A/co < 5M0 

M„„,,bar = < MpeNi + /fb(Af " A/poNi) 5 < McO < 7.6 

I M Mco > 7.6 Mq 

(2) 

where M is the pre-supernova mass of the star, and /fb 
is the fall-back factor, i.e. the fraction (from to 1) of 
the stellar envelope that falls back. The value of /fb is 
interpolated linearly between A/co = 5 Mq (/fb = 0) and 
Mco = 7.6 M0 (/fb = 1). The regimes of no fall-back 
(A/co < 5 M0), partial fall-back (5 < Mco < 7.6 M0) 
and direct collapse (A/co > 7.6 Mq) are estimated from 



core collapse models of Fryer, Woosley & Hartmann (1999) 
and the analysis of Fryer & Kalogera (2001). 

We also allow for NS formation through ECS (e.g., Pod- 
siadlowski et al. 2004). Following Hurley et al. (2000) we 
use the He core mass at the AGB base to set the limits for 
the formation of various CO cores. If the He core mass is 
smaller than A/cburi the star forms a degenerate CO core, 
and ends up forming a CO WD. If the core is more massive 
than A/cbur2 = 2.25 Mq the star forms a non-degenerate 
CO core with subsequent burning of elements until the for- 
mation of a FeNi core which ultimately collapses to a NS 
or a BH. Stars with cores between A/cburi and A/cbur2 may 
form partially degenerate CO cores. If such a core reaches 
a critical mass {Mco,crit = 1-08 Mq, Hurley et al. 2000), 
it ignites CO off-center and non-explosively burns CO into 
ONe, forming a degenerate ONe core. If in subsequent evo- 
lution the ONe core increases its mass to M^cs = 1-38 Mq 
the core collapses due to electron capture on Mg, and forms 
NS (we will refer to a NS formed in this way as ECS NS, 
as opposed to a regular FeNi core collapse compact object 
formation). The ECS NSs are assumed to have unique 
masses of A/rem,bar = -^ecs- If thc ONc corc mass remains 
below A/ocs the star forms a ONe WD. 

Hurley et al. (2000) suggested A/cburi = 1.66 Mq cor- 
responding to A/zams = 6.5 M0 for Z = 0.02. Later calcu- 
lations with the same, but updated evolutionary code (El- 
dridge & Tout 2004a, b) indicated that ECS may occur for 
higher initial masses (A/zams <; 7.5 Mq). For our standard 
model, we adopt A/cburi = 1-83 Mq (A/^ams = 7.0 Mq), 
and this results in ECS NS formation above M^ums = 
7.6 Mq (for masses A/zams = 7.0 - 7.6 Mq the ONe core 
does not reach M^cs and a ONe WD is formed) . It is noted 
that binary evolution through RLOF, may either decrease 
the initial the mass of the ZAMS star required to form the 
corc of mass A/cburi (due to rejuvenation) or increase it 
(due to mass loss). Therefore, binary evolution effectively 
leads to wider initial progenitor mass for ECS NS forma- 
tion. Also metallicity, and wind mass loss influences the 
ECS NS formation range. 

The remnant masses calculated above express the mass 
of baryons, and for NS/BH remnants we convert them to 
gravitational masses (M^cm). We use the quadratic rela- 
tion 

A/rem,bar " A/rcm = 0.075 M^^ (3) 

for neutron stars (Lattimer & Yahil 1989; see also Timmes 
et al. 1996), while for black holes we simply approximate 
the gravitational mass with 

A/rc,„ = 0.9 A/, 

cm, bar (4) 

The resulting remnant mass spectrum covers a wide range 
of masses and is presented in Figure [TJ 

In summary. Helium WDs form from progenitors of ini- 
tial masses in range A/zams <, 0.8 Mq, CO WDs from 
A/zams = 0.8 — 7 Mq, and ONe WD are formed in range 
A/zams = 7—7.6 Mq. Ncutrou stars formed through ECS 

(A/zams = 7.6 - 8.3 Mq) havC mass of A/rem = 1-26 Mq. 

Regular core-collapse (of FeNi cores) NS are formed with 

mass A/,.em = 1-36 Mq (A/zams = 8.3 - 18 Mq); A/,cm = 

1.86 Mq (A/zams = 18-20 Mq) and for higher initial 
progenitor masses NSs form up to adopted maximum 
NS mass: A/NS,max- In our standard model we adopt 
A/NS,max = 2.5 Mq , and above A/zams ^ 21 Mq BHs are 
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arc formed. It is found that single stars may form BHs up 
to ^ 11 Mq for solar metallicity {Z = 0.02) and ^ 30 Mq 
for lower metallicities {Z = 0.001 — 0.0001), which is con- 
sistent with the current observations of the most massive 
BHs in Galactic X-ray transients. The initial-final mass 
relation described above and presented in Figure [T] holds 
only for single stellar evolution and for a specific metallic- 
ity {Z = 0.02). Effects of binary evolution, in particular 
mass loss/gain in RLOF, may alter the initial- final mass 
relation in two ways. First, a compact object may have a 
different mass; higher if a progenitor (or a compact object 
itself) accreted mass or smaller if a progenitor lost mass 
in RLOF. Second, the initial mass limits for formation of 
a given type of compact object are not sharp, since with 
binary mass loss/gain stars of various initial masses may 
form compact objects of a given mass and type. These 
limits are blurred by binary evolution. 

The mass estimates of neutron stars in relativistic dou- 
ble neutron star binaries point to a NS formation mass 
of - 1.35 Mq (Thorsett & Chakrabarty 1999). It is also 
suggestive that the NS mass in Vela X-1 is high ~ 1.9 Mq. 
Since Vela X-1 is a high-mass X-ray binary, system with 
wind-fed accretion (small mass capture efficiency) and a 
massive star donor (short lifetime), the NS probably has 
not accumulated much mass, and the measured mass is 
close to its formation mass. Our adopted model for NS 
formation masses (Timmes et al. 1996) falls in qualitative 
agreement with these observations. Finally we allow com- 
pact object masses to increase through accretion in binary 
systems. Accretion and mass accumulation onto WDs is 
described in detail in §5.7. For NS we need to adopt a 
maximum NS mass, over which NS collapses to BH. Such a 
collapse may lead to a short-hard Gamma-ray burst event. 
Depending on the preferred equation of state the maxi- 
mum NS mass may vary in a wide range (^2 — 3 Mq), 
and in particular may reach ^ 3 Mq if rotation is included 
(Morrison, Baumgarte & Shapiro 2004). At the moment 
the highest measured NS mass is 2.1 ± 0.2 Mq for a mil- 
lisecond pulsar in PSR J0751-f 1807; a relativistic binary 
with helium white dwarf secondary (Nice et. al. 2005). As 
stated above we adopt AfNS.max = 2.5 Mq for our standard 
model, but we relax this assumption in parameter studies. 



2.3.2. Wind 



loss 



The compilation of stellar wind mass loss rates presented 
in Hurley et al. (2000) has been expanded to include mass 
loss from low- and intermediate-mass main sequence stars. 
We have adopted the formulas of Nieuwenhuijzen & de 
Jager (1990) to calculate the mass loss rates for main se- 
quence stars below ~ 8 Mq. Although the mass loss from 
these stars is not large enough to significantly alter the 
evolution of a mass-losing star, it may play an impor- 
tant role in the formation and evolution of wind-accreting 
close binaries. Even with small mass transfer rates charac- 
teristic for the low- and intermediate-mass main sequence 
stars, the X-ray luminosities for accreting BHs and NSs are 
high enough to be detected in deep Chandra exposures. A 
number of faint point X-ray sources were discovered in the 
Galactic center with deep exposures (Wang, Gotthelf & 
Lang 2002; Muno et al. 2003), some of which may be 
explained in terms of wind-fed close binaries (Pfahl, Rap- 
paport & Podsiadlowski 2002a; Bleach 2002; Willcms & 
Kolb 2003; Belczynski & Taam 2004b). 



2.3.3. Rotational velocities 

A compilation of updated observational data on rota- 
tional velocities is used to initiate the stellar spins on the 
ZAMS. The spin evolution is followed as detailed here for 
single stars and in § 3 for binary components. In order to 
obtain a functional form of the relation of the equatorial 
rotational velocity and stellar mass, we used the compila- 
tion of rotational velocities of Stauffer & Hartmann (1986) 
for stars in open clusters. The difference between cluster 
and field stars is quite small for massive stars (with a max- 
imum difference of 10% for intermediate B-type stars), 
but can be as high as 40% for stars later than F-typc, with 
field stars having systematically lower rotational velocities. 

The mean rotational velocity Vrot was determined from 
the projected velocity (vrotsini) assuming a random dis- 
tribution of angles with sini = 7r/4. We fitted Viot as a 
function of stellar mass, and we obtained the following 
empirical functional form 
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and Mo = 6.35j;2;i (errors are at the 1 a level). We stress 
that this is only an empirical functional form of the equa- 
torial rotational velocity as a function of stellar mass. In 
Fig O we present the observational data from Stauffer & 
Hartmann (1986), together with the best fit function. In 
the bottom panel of this figure we also show the ratio of 
the Stauffer & Hartmann (1986) data and the model. 
The spin angular momentum of a star may be expressed 

as 

Ji,spin = ht^i = hMiRfuji (6) 

where, LUi = Vrot /Ri is the angular rotational velocity and 
the coefficient fc; varies as the star evolves and its inter- 
nal structure changes (e.g., it is 2/5 for a solid sphere and 
2/3 for a spherical shell). Following Hurley et al. (2000) 
we consider two structural components for each star: a 
core and an envelope. The spins of these two components 
may decouple in the course of evolution, although we keep 
them coupled in our standard model calculations. The 
spin angular momentum of a star is then 



<^i,spin — [fci,cnv(-^'-^i ~ Mi^q)R^ /Cj., 



(7) 



We use different values than Hurley at al. (2000) 
for the internal structure coefficient fcj. For 
stars with no clear core-envelope structure (A'i = 
0,1,7,10,11,12,13,14,16,17) we use simple polytropic 
models (e.g., Lai, Rasio k. Shapiro 1993) with n = 1.5 and 
n = 3 for low-mass and high-mass objects, respectively, 
giving 



0.205 
0.075 



Mi 



< 1 Mq 
> 1 Mq 







(8) 



For giants with a clear separation between core and en- 
velope (ATi = 3,4,5,6,9) and additionally for stars in the 
Hertzsprung gap {Ki ~ 2,8) we use detailed models of 
giant envelopes (Hurley et al. 2000) and for the core we 
apply a polytropic model with n = 1.5 to obtain 



0.1, fci. 



0.205 



(9) 
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Conservation of the spin angular momentum of a star 
is used then to determine its rotational velocity. Addi- 
tional angular momentum losses from magnetic braking 
(see § 3.2) are also taken into account. 

2.3.4. Convective/ Radiative envelopes 

Stars with convective and radiative envelopes respond 
differently to various physical processes (e.g., magnetic 
braking, tidal interactions or mass loss). Stars that have a 
significant convective envelope are: low-mass H-rich MS 
stars {Ki = 0, 1) within the mass range of 0.35 Mq — 
-^ms.conv, whcrc Mms.conv IS the maximum mass for a 
MS star to develop a convective envelope; giant-like stars 
{Ki = 3, 5, 6) independent of their mass and evolved low- 
mass Helium stars {Ki = 9) below Mhcconv = 3.0 Mq. For 
stars crossing H-rich Hertzsprung gap [K = 2) and core 
helium burning stars {Ki = 4 we obtained detailed mod- 
els using the code described in Ivanova & Taam (2004) 
to check how far from Hayashi line stars cross the border 
between radiative and convective envelopes; stars cooler 
than log(reff) = 3.73 ± 0.02 have convective envelopes, 
while hotter stars have radiative envelopes. We have also 
examined the models presented by Schaller et al. (1992) 
for both solar metallicity and Z = 0.001 and have found 
that the above temperature cut works rather well for these 
two classes of stars for both metallicities. MS stars with 
masses below 0.35 are fully convective. The value of 
M, 



ins,conv 



) depends strongly on metallicity 



M, 



ms.conv 




55.73Z + 0.747 



Z > 0.02 

0.001 <Z < 0.02 

Z < 0.001 



(10) 

Values of Mms.conv in metallicity range Z = 0.001 — 0.02 
are obtained from a fit to detailed evolutionary calcula- 
tions (Ivanova 2006). All other stars (e.g., K = 8 ~ Helium 
stars in the Hertzsprung gap) are assumed to have radia- 
tive envelopes. We use the original Hurley et al. (2002) 
formulas to calculate the mass and depth of convective 
envelopes. 

2.3.5. Helium star evolution 

We assume that low-mass evolved Helium stars {Ki = 9) 
below Mho,conv = 3.0 Mq (as opposed to 2.2 Mq in Hurley 
et al. 2000) expand and form deep convective envelopes in 
their late stages of evolution (e.g., Ivanova et al. 2003; 
Dewi & Pols 2003). Helium stars with convective en- 
velopes are subject to strong tidal interactions (convec- 
tive tides as opposed to radiative damping, see §3.3), 
and if found in an interacting binary, they may alter sig- 
nificantly the fate of a given system. All helium stars 
(Ki = 7, 8, 9) may be subject to stable RLOF. However, 
in dynamically unstable cases we assume a binary compo- 
nent merger in the case of a HeMS donor {Ki = 7) and a 
Helium Hertzsprung gap donor (A'j = 8; e.g., Ivanova et 
al. 2003) or we follow a given system through a CE phase 
for evolved He star donors {K[ = 9 and test whether the 
system survives or merges. The examination of RLOF 
stability and development of dynamical instability are de- 
scribed in detail in § 5). 

The treatment of helium stars is important, for exam- 
ple, in later stages of evolution leading to double neutron 



star formation. The immediate consequences, leading to 
the formation of a new class of close double neutron stars, 
were discussed in Belczynski & Kalogera (2001) and Bel- 
czynski, Bulik & Kalogera (2002a). The formation of the 
new class of (ultracompact) double neutron stars involves 
Helium Hertzsprung gap stars initiating CE phase. Now, 
guided by the better understanding of CE phase survival 
(Ivanova et al. 2003), we do not allow survival in such 
cases, and in our reference model both double neutron 
star merger rates and double neutron star properties have 
changed (Belczynski et al. 2007). However, we note that 
it is still predicted that a significant fraction (~ 40%) of 
close double neutron star binaries form in very close or- 
bits (merger times smaller than 100 Myr). Due to signifi- 
cant updates of the code and new observational results on 
short GRBs with double neutron stars suggested as their 
progenitors (e.g.. Fox et al. 2005) new StarTrack calcu- 
lations relevant to the double neutron star formation have 
been performed (Belczynski et al. 2007; Belczynski et al., 
in prep.). 

3. BINARY ORBITAL EVOLUTION 

Throughout the course of binary evolution we track the 
changes in orbital properties. A number of physical pro- 
cesses may be responsible for these changes. In the general 
case of eccentric orbits we numerically integrate a set of 
four differential equations describing the evolution of or- 
bital separation, eccentricity and component spins, which 
depend on tidal interactions as well as angular momentum 
losses associated with magnetic braking, gravitational ra- 
diation and stellar wind mass losses. For circular orbits 
with synchronized components, we can obtain an exact 
solution for the change of orbital separation using conser- 
vation of angular momentum. Losses of angular momen- 
tum and/or mass associated with RLOF events, magnetic 
braking and gravitational radiation arc taken into account. 

We assume that any system entering RLOF becomes cir- 
cularized and synchronized (if it had not already reached 
this equilibrium state before RLOF). In such a case we 
circularize to the periastron distance 



Ofin = aint(l - e) 
efin = 



(11) 



where and int, fin mark initial and final values), and both 
components are synchronized to the new mean angular or- 
bital velocity. Non circular orbits and non synchronous 
RLOF cases are important for massive binaries, in which 
massive stars may remain eccentric/unsynchronized as 
tidal interactions are not as effective as for low mass stars 
(see §3.3). In particular, the vast majority of close (co- 
alescing) NS-NS binary progenitors may evolve through 
such configuration more than once. It should be noted 
that our prescription (eq. Ilip does not conserve total an- 
gular momentum (i.e., sum of orbital and spin angular 
momenta of the components) . 

For most of the cases in double compact object pro- 
genitor evolution, our prescription leads to a moderate 
(~ 20%) loss of angular momentum. For systems which 
have not been circularized and synchronized before enter- 
ing RLOF there might be substantial mass loss (e.g.. Hut 
& Paczynski 1984), and although this is not taken into ac- 
count in our calculations, it may also lead to some angular 
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momentum loss. 

For some systems, consisting of low mass object (e.g., 
a neutron star) and a massive star (e.g., a massive CHeB 
star) eccentricity is induced by the supernova explosion 
and once the massive star overfills its Roche lobe we cir- 
cularize/synchronize the system, and this may lead to a 
slight increase of total angular momentum (~ 10%). How- 
ever, it needs to be stressed that there is no solution that 
conserves total angular momentum in such cases, as sys- 
tems evolve toward a Darwinian unstable state (see § 3.3) 
and eventually evolve towards a common envelope phase. 
Since the orbital shrinkage of a system with a massive star 
donor during the common envelope phase is usually dra- 
matic (~ 2 orders of magnitude), the moderate change of 
orbital separation prior CE phase (eq. [TT]) and its specific 
magnitude (in the prescription used) does not play an im- 
portant role. 

Violent processes like SN explosions or common en- 
velope phases are taken into account in binary orbital 
evolution. Also nuclear evolution of components (expan- 
sion/contraction affecting stellar spins) is considered. In 
what follows we describe the elements used to calculate 
the orbital evolution. The orbital angular momentum of 
the binary and its mean angular velocity are expressed as 



MiAh^ aG{Ah + Ah) 

Ml + il/2 



(12) 



(13) 



where G is the gravitational constant. 

3.1. Gravitational radiation 

Binary angular momentum loss due to gravitational ra- 
diation is estimated for any type of binary following Peters 
(1964) 



dJgr/ dt 



32 GiMi^MzVAfi + M2 7 



c^a2 (1 — e^)2 



(l + -e2) (14) 



where, c is the speed of light. Emission of gravitational 
radiation causes orbital decay as well as circularization, 
both taken into account during the evolution of a binary 
system. For any given system, the merger time may be 
easily estimated (e.g., see eq.l4 in BKB02). 

3.2. Magnetic Braking 

Each binary component may be subject to magnetic 
braking, causing the decrease of the component's rotation. 
In the case of a detached binary configuration magnetic 
braking is applied directly to the component spins, while 
during RLOF the effects of magnetic braking are applied 
to the orbit, since the components are then kept in syn- 
chronism. Three different prescriptions for magnetic brak- 
ing are incorporated within the StarTrack code and may 
be used interchangeably for parameter studies. In what 
follows we provide a detailed description of the specific 
braking laws adopted. 

Magnetic braking is applied to stars with a significant 
convective envelope, i.e., for low-mass H-ricli MS stars. 



H-rich giant-like stars and cool HG and CHeB stars (see 
§2.3.4 for details) with the exception of low-mass evolved 
Helium stars for which there is not much known about 
magnetic fields. For fully convective MS stars {Ki = 0, 
M < 0.35 M0) magnetic braking may also operate, al- 
though it has been hypothesized that the braking is sup- 
pressed (Rappaport, Verbunt & Joss 1983; Zangrilli, Tout 
& Bianchini 1997) in order to provide an explanation of 
the observed period gap for cataclysmic variables. There- 
fore we assume that magnetic braking is not operative 
for fully convective stars, independent of the prescrip- 
tion used. Since massive core helium burning stars, more 
massive H-rich MS stars, and He-rich MS stars have ra- 
diative envelopes, we assume that magnetic braking does 
not operate in these stars. The prescription for the loss 
of angular momentum associated with magnetic braking 
dJi^mh/ dt takes several forms. Historically, most stud- 
ies have adopted the form suggested by Rappaport et al. 
(1983) where 



dJi^mhl dt = -5.8 X IQ-'^^A'URfuj, 



(15) 



with parameter 7 = 2 in our model calculations. However, 
studies based on the observations of rapidly rotating stars 
show that the Skumanich relation (J ex w^) is inadequate 
in this regime and point to a weakening of magnetic brak- 
ing due to saturation of the dynamo (Andronov, Pinson- 
neault & Sills 2003). In this case, the angular momentum 
loss rate takes the form 



dJi^mb/dt = -8.88x10^^2 V^R^T^ 



LOi UJi < Wcrit 

2 

'^■i'^crit 'J-'crit 

(16) 

where, i denotes the component for which magnetic brak- 
ing is operating, uji [ Myr~^] is angular velocity, and cjcmt 
stands for a critical value of angular velocity above which 
the angular momentum loss rate enters the saturated 
regime. If the latter law is used, the saturation is applied 
only for MS stars and cJcrit is interpolated from Tabic 1 of 
Andronov et al. (2003). 

In addition, we also include the form of magnetic brak- 
ing from the results of a study by Ivanova & Taam (2003). 
In this latter study, an intermediate form of the angular 
momentum loss rate was derived (J oc u^'^) based on a two 
component coronal model as applied to the observational 
data relating stellar activity to stellar rotation. Specifi- 
cally, we adopt 



dJi^mh/ dt = -619.2i?j^ 



(w.,/9.45 X IfS^f 



10i-^(wi/9.45 X lO^y-^ UJi > uj^ 

(17) 

with Wx = 9.45 X 10® Myr ^. This law is used for the 
StarTrack standard model calculations. 

3.3. Tidal Evolution 

The evolution of the orbital parameters (a, e) as well 
as component spins (wi, i = 1, 2) driven by tidal inter- 
actions of binary components is computed in the stan- 
dard equilibrium-tide, weak- friction approximation (Zalin 
1977, 1989), following the formahsm of Hut (19810- This 
formalism allows us to treat binaries with arbitrarily large 



^Note that upon entering RLOF any binary system is instantly synciironized and circularized. 
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eccentricities. We assume that the only sources of dissi- 
pation are eddy viscosity in convective envelopes and ra- 
diative damping in radiative envelopes. Specifically, we 
integrate numerically the following differential equations 
in parallel with the stellar evolution 



with i?i.cnv denoting the depth of the convective envelope 
and Li the bolometric luminosity of a given component 
(Rasio et al. 1996). 

The numerical factor /i^conv is defined as 



fijCoiYV — TTlifl 



1, 



Pi 



tid 



2ri,, 



^j^^g2'j _ ^-^ _ ^2-^3/2 y^(c^) ^ (18) with the tidal pumping timescale Pi^tid defined as 



^) = -27^^tid(|),*(l + gO(f)'TTr^ 



tid 



/3(e2)-ii(l-e2)3/2^4(e2):f- (19) 



du>i 
~dt 



tid 



3Ftid (t); ^ 



(Ik\'^ "orb 

I a J (l-e^)« 



(/2(e2)-(l-e2)3/2/,(e2)^) (20) 



where 



/2(e2) 
Me') = 1 
/4(e2) = 1 
Me' ) - 1 



2 " ' 8 ' 16 

l + fe' + f + AgS 



15 2 
4 ^ 

2"^ ^ 8 
.2 _L 3g4 



15 „4 
8 

o ^ 



Ap6 

64 



3e2 



and Ti^gyr is the gyration radius and is defined by /j = 
Mi(ri_gyri?i)^, with li denoting the moment of inertia of a 
given binary component. Here the mass ratio is defined as 
follows, 

_ / M2/M1 I = 1 
M1/M2 i = 2 



9i 



(21) 



The quantity {k/T)i is the ratio of the apsidal motion 
constant k (which depends on the interior structure of the 
star) over the timescale T of tidal dissipation. Following 
Hurley, Tout & Pols (2002), we calculate that constant 
for either the equilibrium tide with convective damping 
{{k/T)i = (fc/r)i.con) or the dynamical tide with radia- 
tive damping {{k/T)i = (A:/T)i_rad)- Radiative damping is 
applied to stars with radiative envelopes: MS stars with 



mass above My, 



CHeB stars with mass above 7 Mq 



massive evolved He stars and all He MS stars. For all other 
stars, convective damping is applied (see § 2.3.4 for details 
on convective/radiative envelopes). We do not calculate 
tides on stellar remnants, e.g., on WDs {Ki > 10). 

The constant for convective damping is obtained from 



2 M 



Mi, 



21 Ti, 



Mi 



cnv — 1 



(22) 



where Mj^env is the mass contained in the convective en- 
velope of component i. The eddy turnover time Ti conv is 
calculated as 



0.431 



-^^i , en V -^i , cnv ( -^i o : ) 



3Li 



1/3 



i,tid 



P 



orb 



P 



i,spin 



(24) 



(25) 



where Porb and Pi. spin sltc the binary orbital period and 
the spin period of component i, respectively. This factor 
represents the reduction in the effectiveness of eddy vis- 
cosity when the forcing period is less than the turnover 
period of the largest eddies (Goldreich & Keeley 1977) 
The constant for radiative damping is calculated from 



rad 



/ M- R -^ 

1.9782 xl0S/^^(l + gi)5/6i;2 yr~' (26) 



where a second-order tidal coefficient E2 = 1.592 x 
lO-^Mj^ **"^ was fitted (Hurley et al. 200^1) to values given 
by Zahn (1975). 

Finally, we have introduced an additional scaling factor 
Ptid in the evolution equations (eq. [TSl [111 HOI) which we 
normally set to: 



Ptid 



Pt 



tid, con 
id, rad 



50 convective cnv. 
1 radiative env. 



(27) 



and the distinction between the stars with convective and 
radiative envelopes is given in § 2.3.4. This factor makes 
tidal forces (both in case of convective and radiative damp- 
ing) more effective than predicted by the standard Zahn 
theory. The choice of this specific value of Ptid is a result 
of our calibration against the cutoff period for circulariza- 
tion of binaries in M67 and from the orbital decay of the 
high mass X-ray binary LMC X-4 (for details see § 8.2). 

The orbital angular momentum change associated with 
tides is calculated from 



Note that there is a typo in their eq.(42); missing ^. However, 
communication. 



dJiXidl dt — 3Ptid/i 

x(/2(e2)-(l-e2)3/2/5(e2)^) (28) 
and the change of binary parameters is calculated with 

eqs. m [m Eol 

Pre-main sequence tidal synchronization and circular- 
ization. We also allow for pre-MS tidal interactions. Since 
we do not follow pre-MS evolution, all binaries with or- 
bital periods shorter than 4.3 d (Mathieu et al. 1992) are 
simply assumed to have circularized and all binary com- 
ponents to have synchronized by the time they reach the 
ZAMS. For binaries with longer orbital periods we apply 
our assumed distribution of initial eccentricities (see § 7) 
, , and initial rotational velocities for binary components (see 
yr (23) §2.3.3). 

their binary code utiUzes the proper formula. Jarrod Hurley, private 
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Darwin instability. One important consequence of tidal 
interactions in massive binaries is the possible occurrence 
of the Darwin instability (e.g., Lai et al. 1993). When 
the more massive component is spinning slowly compared 
to the orbital rate of its companion, tidal forces will tend 
to spin it up, leading to loss of orbital angular momen- 
tum (orbital decay). Usually this orbital decay will stop 
when synchronization is established. However, if, in the 
synchronized state, more than a third of the total binary 
angular momentum would be in the component spins, then 
synchronization can never be reached and the components 
will continue to spiral in. We follow this process until one 
of the binary components overflows its Roche lobe. In all 
cases RLOF is found to be dynamically unstable (§ 5.1 and 
§ 5.2) and the system goes through a CE phase leading ei- 
ther to a merger, or further orbital decay with envelope 
ejection (§ 5.4). 

3.4. Mass and Angular Momentum Loss from Binaries 

Mass lost from the binary components in stellar winds 
carries angular momentum, in turn affecting the orbit 
through tidal coupling. Similarly, during RLOF, some of 
the transferred material and its associated angular mo- 
mentum may be lost from the system. In this section we 
consider the amount of angular momentum loss associated 
both with stellar winds and RLOF phases. However, for 
RLOF we only consider here dynamically stable phases, 
while the change of the orbit following unstable RLOF 
(common envelope events) is described in § 5.4. 

For stellar winds we assume spherically symmetric mass 
loss, which carries away the specific angular momentum of 
the mass-losing component (Jeans mode mass loss). The 
corresponding change of the orbit (Jeans-mode mass loss) 
is calculated from 



a(Mi -I- M2) 



const. 



(29) 



The above approach holds for circular orbits, however the 
change in binary separation a is similar for eccentric or- 
bits (Vanbeveren, Van Rensbergen & De Loore 1998, see 
p. 124). 

In the case of stable RLOF with compact accretors 
(WD, NS, BH; A'acc = 10,11,12,13,14,16,17) we limit 
(although this assumption may be relaxed) accretion to 
the Eddington critical rate 



Mcdd = 2.088 X 10 



-3 Race 



e{l + X) 



Mo yr-i (30) 



and the corresponding critical Eddington luminosity may 
be expressed as 



GAfaccA^cdd 
Race 



(31) 



where i?acc denotes the radius at which the accretion onto 
compact object takes place (a NS or a WD radius, and 
three Schwarzschild radii for a BH), X denotes the com- 
position of accreted material (0.7 for the H-rich material, 
and 0.0 for all other compositions), and e gives the conver- 
sion efficiency of gravitational binding energy to radiation 
associated with accretion onto a WD/NS (surface accre- 
tion e = 1.0) and onto a BH (disk accretion e = 0.5). We 



also note that above some critical (very high) accretion 
rate, nuclear burning will start on the WD surface. This 
can be much more radiatively efficient than the gravita- 
tional energy release and the above relations break down. 
If the mass transfer rate is higher than Afodd we expect 
the excess material to leave the system from the vicinity 
of the accreting object and thus to carry away the specific 
angular momentum of the accretor. The angular momen- 
tum loss associated with a given systemic mass loss in a 
RLOF phase is obtained from 



dJRLOp/ dt — RI 



rb(l-ia)Md 



(32) 



where i?com — aAldon/(A-/^don + A^acc) is the distance be- 
tween the accretor and binary center of mass, and Afdon 
is the mass transfer rate (donor RLOF rate, see eq. |49|) . 
The /a fraction of material transferred from the donor 
is accreted on the compact object. If mass transfer is 
sub-Eddington then fa = 1 (conservative), otherwise it 
is /a = Medd/Afdon (non-couservative evolution). Here 
we assume that the radiative efficiency is not a function 
of the mass transfer rate. Some work has suggested that 
at high transfer rates, flows onto black holes may become 
radiatively inefficient as photons are trapped in the flow 
and advccted into the black hole (see e.g., Abramowicz 
et al. 1988), or that substantially super-Eddington accre- 
tion may be possible in non-spherical accretion flows (e.g., 
Begelman 2002). In the current version of the code, we do 
not consider these possibilities. 

For all other, non-degenerate accretors {Kacc = 
0, 1, 2, 3, 4, 5, 6, 7, 8, 9) we assume a non-conservative evo- 
lution through stable RLOF, with part of the mass lost 
by the donor accreted onto the companion (/a), and the 
rest (1 — /a) leaving the system with a specific angular 
momentum jioss expressed in units of 27ra^/Poib (see Pod- 
siadlowski. Joss & Hsu 1992). The angular momentum 
loss is then estimated from 



RLOF 



/ dt = jio 



J, 



orb 



"A/do 



A/a„ 



-(l-.fa)Afdo 



(33) 



For our standard model calculation we adopt jioss = 1, 
and /a = 0.5 (half of the transferred mass lost from sys- 
tem, e.g. Meurs & van den Heuvel 1989). However, we 
note that the amount of mass loss (as well as the specific 
angular momentum with which the mass is lost from the 
binary) may change from case to case. Ideally, one would 
like to know the details of mass transfer/loss for all po- 
tential binary configurations, and change jioss and /a ac- 
cording to the types of interacting stars as well as binary 
properties. Since, such predictions or understanding are 
not available, we treat jioss and /a as parameters, which 
are applied evenly to all the stars in a given simulation. 

4. WIND MASS LOSS/ACCRETION IMPLEMENTATION 

We adopt the compilation of mass loss rates from Hur- 
ley at al. (2000). We have further extended the original 
formulas to include winds from low- and intermediate- 
mass MS stars. The structure of the star (and its sub- 
sequent evolution) in response to stellar wind mass loss 
is self-consistently taken into account with the Hurley et 
al. (2000) evolutionary formulas. The most important ef- 
fects include possible removal of the H-rich envelope of a 
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massive star or a more gradual nuclear evolution with de- 
creasing mass. The effects of wind mass loss from binary 
components on the orbital parameters are also accounted 
for (sec § 3.4). 

The effects of mass increase of binary components due 
to accretion from the companion winds are neglected. Ei- 
ther the wind accretion rates are very low or the high 
wind accretion phases do not last for long, which does 
not translate into significant mass increase of a compan- 
ion star. However, we estimate the wind accretion rates 
onto NSs and BHs since it may give rise to bright X-ray 
emission (see § 9) . 

The wind accretion rate is calculated in the general case 
of eccentric orbits, i.e. we obtain accretion rate (and accre- 
tion luminosity) for a specified position on the orbit, or we 
integrate over a specific part of the orbit (e.g., correspond- 
ing to the exposure time of given observations). This may 
be of importance for eccentric wind-fed binaries, e.g., high 
mass X-ray binaries (see § 9.2). We have also implemented 
an orbital-averaged solutions. The two solutions may be 
adopted as required for a given project or analysis. 

4.1. General Eccentric Orbit Case 

We follow the Bondi & Hoyle (1944) accretion model 
to calculate the accretion from stellar wind. As an ap- 
proximation we may express (Boffin & Jorissen 1988) the 
accretion rate as 



acc,wind — Q^wind 



27r(GA /acc)' 



-13/2' 



(34) 



where awind = 1.5 is the accretion efficiency in the Bondi- 
Hoyle model, although it may be as low as 0.05 in some 
specific cases (e.g., see hydrodynamical simulations of Tlie- 
uns. Boffin & Jorissen 1996 for Barium star formation), 
Cwind is the wind sound speed, and K-ei is the relative 
velocity of the wind with respect to the accreting star. 
The local (undisturbed) density of the wind matter p in 
the vicinity of the accreting object may be calculated in a 
steady spherically symmetric case from 



don, will! 



d = -47rHpy, 



wind 



(35) 



where Mdon.wind is the wind mass loss rate form the donor, 
r is the instantaneous distance between the two stars, and 
Vwind is the wind velocity. We assume that the wind fiow 
is supersonic (V^ci ^ Cwind) so that c^i„d may be dropped 
from eq. 1341 We introduce p (expressed through eq. 
into eq. [M] to obtain 



M; 



(GMacc)' 



acc.wmd 



Q^wind (-)T/3 T/ 2 '^^fiori^wind 

rcl ^wind^ 



(36) 



The accretion rate calculated with eq. [32] varies as the ac- 
creting object moves in its orbit around the mass-losing 
star. The relative distance r of the two stars is obtained 
through the Kepler equation for a given orbit. Obviously r 
is a function of orbital position. The vector of the relative 
velocity V^ei is defined as 



^■cl = V'accorb + K 



wind 



(37) 



where I4cc,orb denotes the instantaneous velocity of the 
accretor on the orbit relative to the mass losing star, and 



is readily obtained for a given position through the Ke- 
pler equation. The direction of the wind velocity vector 
Vwind follows tlic vcctor pointing toward the accretor on 
its relative orbit around the mass-losing star. We set the 
wind velocity proportional to the escape velocity from the 
surface of the mass-losing star 



wind 



2/3, 



GM, 



don 



'wind " 



'don 



(38) 



and vary /3wind with the spectral type of the mass- 
losing star. For extended (i?don > 900 R©) H-rich gi- 
ants {Kdon = 2, 3, 4, 5, 6) slow winds are assumed /3wind = 

0. 125. For the most massive MS stars (> 120 Mq) /3wind = 
7, for low mass MS stars (< 1.4 Mq) /3wind = 0.5 and 
the value of /3wind is interpolated in-between. For He-rich 
stars (i^don = 7,8,9); /3wind = 7 for Mdon > 120 M©, 
/3wind = 0.125 for Mdon < 10 Mq, and is interpolated in- 
between. The values of /3wind follow from the observations 
of wind velocities for different type of stars (Lamers, Snow 
& Lindholm 1995; Kucinskas 1999) and are adopted from 
the discussion of wind properties in Hurley et al. (2002). 

4.2. Orbit-averaged Case 

We use eq. [35] to obtain the orbit-averaged accretion 
rate. The wind velocity vector is assumed to be perpen- 
dicular to the orbital speed vector (as on a circular orbit), 

1. e., = Kccorb + Vwind- The wind velocity is taken 
from eq. [35] The orbital velocity of the accretor is taken 
to be constant and is obtained from the circular orbit ap- 
proximation l/^^^_orb = G(Afacc + Mdon)/a- Finally, l/r^ is 
substituted in eq. [35] with its mean value over one orbital 
revolution, i.e., l/(a^\/l — e^) to obtain 



M, 



acc.wmd 



2 

wind 



•^wind -^-^don.wind 



2a2 



(1 



^2)3/2 

(39) 

where Fwind is a parameter (see below) and V"^ = 

acc,orb/ wind 

For highly eccentric orbits, the averaged (over one orbit) 
accretion rate calculated with the eq. [35] may exceed the 
companion mass loss rate. This is a direct result of the or- 
bital averaging used above. To avoid this we follow Hurley 
et al. (2002; §2.1) and adopt Fwind such that Mace, wind 
never exceeds 0.8Mdon,wind- 

5. ROCHE LOBE OVERFLOW CALCULATIONS 

Different physical processes may be responsible for driv- 
ing RLOF. In the following we describe the treatment of 
mass loss and mass accretion in our model. 

5.1. Mass Transfer/ Accretion Rate 

For any binary system during RLOF phases with a non- 
degenerate donor (ifdon < 10) we calculate the radius mass 
exponents for the donor and its Roche lobe 



don 



Uob 



_ d In J?don 

d In Afdon 

d In -RdonJob 

d In Mdon 



(40) 



(41) 



10 



and wc estimate the change of donor radius with time due 
to its nuclear evolution as 



Covl — 



d In Rdo 
dt 



(42) 



The above derivatives are calculated numerically with the 
use of the analytic single star formulas of Hurley et al. 
(2000). The time derivative of the stellar radius (Cevi) is 
obtained directly from the single star formulas, since the 
radius of a given star is tracked in time. To obtain re- 
sponse of the donor to mass loss (Cdon) we calculate the 
response of the star to an instantaneous (over a timestep 
of only 1 yr, a time interval unimportant for stellar evolu- 
tion) mass loss through (artificially) increased wind mass 
loss. Finally, the Roche lobe exponent is obtained by re- 
moving 1% of the donor mass, part of which is transferred 
to the accretor and the rest is lost with the specific angu- 
lar momentum of the accretor from the binary (see §3.4); 
a new Roche lobe radius and the numerical derivative are 
then readily calculated. 

RLOF may be driven by different physical processes; 
angular momentum losses connected to magnetic brak- 
ing (applied directly to the orbit, given the assumption 
of synchronism during RLOF) and gravitational radiation 
or expansion due to nuclear evolution. Wc do not include 
tides as we assume that binary is circular and synchro- 
nized during RLOF. The timescales for magnetic braking, 
and gravitational radiation are calculated from 



Tmh 



orb 



^^don.nib/ dt -|- dJsLCcmh / dt 
'J orb 



igr 



djgi/ dt 



(43) 



(44) 



where expressions for dJg^/ dt, dJimb/ dt, are given in 
§3.1, and §3.2 respectively. 

If RLOF is driven by the combination of angular mo- 
mentum losses changing the orbit and nuclear evolution of 
the donor we then calculate the mass transfer rate from 



Covl + :;r- + :f- 

Meq = ^ ""^ "' j\/don 



Cdon Clob 

and the corresponding mass transfer timescale 



(45) 



(46) 



Additionally wc estimate the thermal timescale for the 
donor following Kalogera & Webbink (1996) from 



30 X Maon 



-^don-^don 

and the mass transfer rate on the thermal timescale 



Mth = — 



(47) 



(48) 



In the case of stable RLOF, Tcq > rth, a donor is in ther- 
mal equilibrium, and we use eq. to calculate the mass 
transfer rate. Otherwise, for r^q < Tth, RLOF proceeds 
on the thermal timescale and we evolve a given system 



calculating the mass transfer rate from eq. HH Wc follow 
the timescales of the donor as it evolves through RLOF, 
and apply the appropriate mass loss rate. For example, 
a massive donor may be transferring mass on the thermal 
timescale at first, but once it loses a fraction of its mass, 
the mass transfer becomes stable and RLOF proceeds on 
the timescale defined by Tcq. However, in some cases the 
RLOF is so rapid that it may eventually lead to a dy- 
namical instability. Once M^q changes sign and becomes 
positive, the donor loses its equilibrium, and the system 
evolves either on the thermal or dynamical timescale. In 
this case a special diagnostic diagram is used (see below) to 
decide which of the two timescales is relevant. We also al- 
low for the development of a delayed dynamical instability, 
which may occur for stars with a radiative envelope, but 
with a deep convcctive layer. Dynamical instability during 
RLOF leads to a spiral-in of the binary components and 
common envelope evolution (CE). Wc follow the CE phase 
to determine whether the binary survives (ejection of the 
envelope at the expense of orbital energy) or if a merger 
of the binary components (single star formation) occurs. 

The following summarizes the calculation of the RLOF 
mass transfer rates 

CE /merger Mdon > gddi x Mace 

Afoq Meq < and Tcq > Tth 

Mth Meq < and Tcq <= Tth 

Mth/ C E /merger diagnostic diagram 

(49) 

where we additionally assume that above some critical 
mass ratio (gddi = -A/don /-A/acc) the binary system will 
evolve toward delayed dynamical instability (Hjellming & 
Webbink 1987), leading to rapid inspiral and CE evolution. 
For H-rich stars Hjellming (1989) gives a range q^di = 2 — 4 
depending on the evolutionary state of a donor, while 
Ivanova & Taam (2004) obtain q^di = 2.9 — 3.1. In our 
standard model calculations we adopt qddi = 3 for H-rich 
stars (A'i = 0,1,2,3,4,5,6). For He-rich stars we adopt 
critical mass ratios from Ivanova et al. (2003); qddi = 1-7 
for HcMS stars {Ki = 7), while ^ddi = 3.5 for evolved He 
stars {Ki = 8, 9). We note that the study of Ivanova et al. 
(2003) was targeted for He stars with NS accretors only. 
However, we adopt their values for systems with He star 
donors and arbitrary accretors, since detailed models for 
arbitrary accretors are not available. Also, dynamical in- 
stability may be encountered if the trapping radius of the 
accretion flow exceeds the Roche lobe radius of the accre- 
tor (§5.4). Additionally, we consider the case of spiraMn 
in the case of the Darwin instability, where the compo- 
nents' spin angular momentum is comparable to the or- 
bital angular momentum (§ 3.3). 

For the donor stars without a well defined core-envelope 
structure (A'don = 0, 1, 7, 10, 11, 12, 16, 17) we assume that 
dynamical instability during RLOF always leads to a 
merger. The same is assumed for the donors in the 
Hertzsprung gap (i^don — 2, 8) as there is no clear entropy 
jump at the core-envelope transition (Ivanova & Taam 
2004; Belczynski et al. 2007). In the case of a merger a sin- 
gle stellar object is formed. However, we do not follow its 
evolution here, as the chemical composition and structure 
of merged remnants is not well understood and certainly is 
different than normal stars. This may lead to an underes- 
timate of our synthetic supernovae rate, since potentially 
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some merger products are massive enough to evolve and 
explode as Type II or Ib/c SNe. For H-rich and He-rich 
giant-like donors (ifdon = 3, 4, 5, 6, 9) we follow CE evolu- 
tion, and assuming ejection of the entire donor envelope, 
we calculate the most probable outcome with conservation 
of energy (see § 5.4). If RLOF is encountered for a system 
with an evolved Helium star donor (A'i = 8, 9), then it is 
found that for low donor masses (< 4 — 5 Mq) RLOF is 
stable (although it may proceed at very high rates) while 
for higher donor masses it leads to a CE phase (e.g., see 
Ivanova et al. 2003). The survival of the binary then de- 
pends on the donor properties (e.g., stellar structure, en- 
velope binding energy, its mass, binary separation) and in 
particular for He stars in the Hertzsprung gap (i^don = 8) 
CE phase always leads to merger. 

The mass accretion rate in a dynamically stable RLOF 
is calculated from 

Mace = fj'idon (50) 

where Mdon is the donor RLOF mass transfer rate (see 
eq. The parameter denotes the fraction of the 

transferred mass which is accreted, while the rest (1 — /^) 
is ejected from the system (see §3.4). Mass accretion in 
dynamically unstable cases (CE events) is calculated only 
for NS and BH accretors, since only then significant accre- 
tor mass gain may be expected in spite of the very short 
timescales (for details see BKB02). 

5.2. Diagnostic Diagram for Rapid Mass Transfer 

The aforementioned diagnostic diagram is shown in 
Fig [31 Once RLOF proceeds on the thermal timescale, 
and the donor is no longer in thermal equilibrium, we do 
not have proper stellar models to use and calculate the 
donor properties (e.g., RLOF rate). Therefore, we use an 
approximate method and calibrate it based on the results 
from detailed stellar evolutionary and mass transfer cal- 
culations, which arc not limited to stars in thermal equi- 
librium. When the donor loses its equilibrium, we use the 
stellar and binary properties to predict whether the sys- 
tem will evolve through the phase of thermal mass transfer 
and the donor will regain its equilibrium, or the RLOF will 
become dynamically unstable and will eventually lead to 
CE evolution. We plot the donor Roche lobe radius versus 
decreasing donor mass under the assumption that mass 
transfer is non-conservative and proceeds on the thermal 
timescale (see eas. l47land |48)) . For NS/BH accretors the 
accretion rate is limited by the Eddington rate, while for 
all other accretors, a fraction /a of transferred material is 
accreted. The associated specific angular momentum loss 
is described in §3.4. As the mass of the donor decreases 
with mass transfer the Roche lobe first shrinks and then at 
some critical mass ratio (giow), it starts expanding again 
(see the solid line on the top panel. Fig. 2). If the mass 
ratio at the moment the star loses its equilibrium qi^t is 
not greatly different than giow we expect that the donor 
may regain the equilibrium when the system is expanding. 
The dashed line arrow in Figure 2 shows the expected be- 
havior of the donor when it loses its equilibrium. If the 
system does not evolve into a CE phase then we expect the 
donor to regain its equilibrium at the position indicated by 
the arrow. Of course this is just an approximation, since, 
as the donor evolves, the radius-mass exponent changes. 



We use a number of published (Tauris & Savonije 1999; 
Wellstein & Langer 1999; Wellstein, Langer & Braun 2001; 
Dewi & Pols 2003) and unpublished (N. Ivanova 2004, pri- 
vate communication) detailed calculations to calibrate the 
diagnostic diagram. Based on these studies we find that a 
CE phase ensues if 

/ 9mt > 1-2 giow -f^don = 2, 3, 4, 5, 6 
\ gint > 2.0 giow IQon = 0, 1, 7, 8, 9 

Otherwise the system is evolved through RLOF on the 
donor's thermal timescale. 

5.3. Thermal Timescale Mass Transfer 

Once a binary is identified as a thermal timescale RLOF 
system, we assume that the mass transfer rate remains 
constant throughout the entire episode. We calculate the 
rate using eq. 25] where we use properties corresponding 
to the time the donor loses its thermal equilibrium. This 
may be justified by the following: (i) thermal mass trans- 
fer rates have been shown to be rather constant within 
a factor of ^ 2 — 3 (Paczynski 1971), (ii) since the rates 
are calculated at the time the star loses equilibrium, it is a 
good approximation (and the best possible with only equi- 
librium stellar models being available) for the short lived 
phase of thermal mass transfer that follows. 

In the bottom panel of Figure [3] we show an example 
calculation through a thermal RLOF phase, followed with 
a slower (driven by nuclear evolution) RLOF period af- 
ter the donor has regained its thermal equilibrium. The 
specific system was chosen to match the RLOF calcula- 
tion of Wellstein et al. (2001) for a 16 Mq and 15 Mq 
binary with an initial period of 8 days. The RLOF 
starts when the primary evolves off the main sequence 
and crosses the Hertzsprung gap. Mass transfer initially 
proceeds on a thermal timescale at a very high rate (~ 
2.8 X 10^"^ Mq yr~^), then the star regains its equilibrium 
and the RLOF rate decreases with time by more than or- 
der of magnitude (~ 10~* Mq yr~^). Our calculation can 
be directly compared to Wellstein et al. (2001): see their 
Figure 4, left panel. Their detailed stellar evolution cal- 
culation shows a thermal RLOF rate of ^ 10^"^ Mq yr~^, 
followed by a slower RLOF phase characterized by rates 
of ~ 10"'' Mq yr~^, very similar to what we find with our 
simplified prescription. Our RLOF phase lasts about twice 
as long as that of Wellstein et al. (2001), who in contrast to 
our calculation assumed conservative evolution and did not 
include effects of tidal spin-orbit interactions. We choose 
not to modify our standard model assumptions (e.g., ne- 
glect tidal interactions) for comparisons, and therefore em- 
phasize some differences with previous calculations. More 
comparisons of RLOF sequences are presented in §8.1. 

5.4. Dynamical Instability and Common Envelopes 

Dynamically unstable mass transfer may be encountered 
in a number of ways. Most often it is the direct conse- 
quence of stellar expansion during rapid nuclear evolution 
phases. However, loss of orbital angular momentum (e.g., 
via magnetic braking, gravitational radiation, or tides) 
may also lead to dynamical instability. 

Additionally, we allow a system to evolve into a CE 
phase if the trapping radius of the accretion flow exceeds 
the Roche lobe radius of the accretor. The trapping radius 
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is defined as (Bcgelman 1979) 



ap 



Mdon i?a 



(52) 



odd 



proposed by Nelemans & Tout (2005), based on the non- 
conservative mass transfer analysis by Paczynski & Zi- 
olkowski (1967), with the assumption that the mass loss 
reduces the angular momentum in a linear way. This leads 
to reduction of the orbital separation 



Following King & Begelman (1999) and Ivanova et al. 
(2003) we check whether the mass transfer rate exceeds 
a critical value above which the system is engulfed in a 
CE 

r» Ti'r -Race. lob /ro\ 

2 X Medd^T- — (53) 



where i?acc.iob is the accrctor Roche lobe radius, and Mcdd 
is the Eddington critical accretion rate (see eg. [50)1 . 

Below we present two different implementations of the 
orbital contraction calculation during CE that arc incor- 
porated in StarTrack. 

Standard Energy Balance Prescription. If dynamical 
instability is encountered a binary may enter a CE phase. 
We use the standard energy equation (Webbink 1984) to 
calculate the outcome of the CE phase 



GAfdon.fin-^'^acc GAfjon.int-^'^a 



2A, 



GM, 



don , int A^don ,cnv 



..lob 

where, A'/don.onv is the mass of the donor envelope ejected 
from the binary, i?don,iob is the Roche lobe radius of the 
donor at the onset of RLOF, and the indices int, fin denote 
the initial and final values, respectively. The parameter A 
is a measure of the central concentration of the donor (de 
Kool 1990; Dewi & Tauris 2000). The right hand side of 
equation (j54p expresses the binding energy of the donor's 
envelope, the left hand side represents the difference be- 
tween the final and initial orbital energy, and ctco is the 
CE efficiency with which orbital energy is used to unbind 
the stellar envelope. If the calculated final binary orbit is 
too small to accommodate the two post-CE binary com- 
ponents then a merger occurs. In our calculations, we 
combine ctcc and A into one CE parameter, and for our 
standard model, we assume that aco x A = 1.0. This is for 
all but evolved naked Helium stars {Ki = 9) for which we 
adopt Qfco = 1.0 and A = O.Si?;^"'®, where Ri is radius of 
Helium star in solar radii. The relation for A was obtained 
with Ivanova's (2003) evolutionary code. If a compact ob- 
ject spirals in the common envelope it may accrete signif- 
icant amounts of material because of hyper-critical accre- 
tion (Blondin 1986; Chevalier 1989, 1993; Brown 1995). 
We have incorporated a numerical scheme to include the 
effects of hyper-critical accretion on NSs and BHs in our 
standard CE prescription (for details see BKB02). Com- 
pact objects gain, on average, several tenths of solar mass 
in CE if hyper-critical accretion is allowed. However, we 
also allow for evolution with no hyper-critical accretion 
following recent results of accretion flow calculations with 
geometry specific for compact object moving through com- 
mon envelope. These calculations indicate that accretion 
can be limited to only 0.01 Mq (E.Ramirez-Ruiz, private 
communication) . 

Alternative Angular Momentum Prescription In addi- 
tion to the standard prescription for common envelope 
evolution based on comparing the binding and orbital en- 
ergies (see above), we investigate the alternative approach 



-^fin / -, -A^don.onv 

= 1-7 



tot, fin 



Aftot.int / A-ftot,int V Mlon,finAfacc,fin 



don, int -^^acc, int 



(55) 

where Mdon.cnv is the mass of the donor envelope lost by 
the system, Mtot,int, A'/tot.fin sue the total masses of the 
system before and after CE, and 7 is a scaling factor. Fol- 
lowing Nelemans & Tout (2005) we use 7 = 1.5 and note 
that hyper-critical accretion is not included in this pre- 
scription. 

The two above prescriptions are extended (e.g., BKB02) 
to the case where both stars lose their envelopes, which 
happens if the stars have giant-like structure {Ki = 
3, 4, 5, 6, 9) at the onset of CE phases (see Bethe & Brown 
1998). 

5.5. Mass Transfer from Degenerate Donors 

Degenerate donors (A'don — 10,11,12,16,17), are also 
considered. The RLOF is assumed to be driven by gravi- 
tational radiation only 



M, 



don 



MdonD 



djgi-/ dt 

Jnrh 



with 
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(56) 



(57) 



where the mass ratio is defined as q ~ A/acc/Afdon, /a de- 
notes the fraction of transferred material that is accreted 
by the companion (defined and evaluated in §3.4), and 

/3mt =M|„„/(A/don+Macc)'. 

5.6. Effects of Mass Transfer on Stellar Evolution 

Mass loss/gain changes the subsequent evolution of 
stars. We implement RLOF mass loss/gain by adding 
an extra term in the original Hurley et al. (2000) stellar 
evolution formulae. In case of mass loss we increase the 
wind mass loss rate to match the combined effects of wind 
and RLOF mass loss. To treat mass gain and potential 
accrctor rejuvenation, we add the RLOF mass accretion 
rate, as calculated in § 5.6, to the accrctor wind mass loss 
rate (they have opposite signs). Additionally, in case of 
rejuvenation evolutionary timescales and stellar ages are 
modified as suggested by Tout et al. (1997) and Hurley et 
al. (2002) by the relative mass change (see the following 
equation) to calculate the net effect on the star. For main 
sequence stars we can calculate the change of the age of a 
given star (due to accretion or mass loss) from 



/r 



■cj ''age, int 

'^ms.int 



(58) 



where is the main sequence lifetime, and indices 
int, fin mark the state before and after some amount 
of mass is transferred, respectively. The factor /roj is 
unity for all mass losing stars and for hydrogen MS stars 
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(A'i = 1) with radiative cores (0.3 < Mi < 1.25 M©), while 
it is /rcj = -Mijnt/-^i,fln for hydrogen MS stars with con- 
vective cores and hchum MS stars {Ki = 7; they have con- 
vective cores) and it reflects the effects of additional fuel in 
the core. For HG stars {K = 2) we change the timescales 
using /rej = 1 as for MS stars with radiative cores. In 
this way we ensure that the subsequent evolution of the 
donor is followed consistently, i.e., evolutionary timescales 
and physical properties of mass losing/gaining stars (e.g., 
core masses) are changed in agreement with stellar mod- 
els. Our wind mass loss formulae are implemented the 
same way as in the original Hurley et al. (2000, see their 
§ 7.1) and the above scheme allows for appropriate changes 
of evolutionary timescales and core masses, both in cases 
of mass loss and gain. 

For simplicity, we assume that the composition of the ac- 
creted material matches that of the accretor, although this 
may not always be the case. Only in the case of accretion 
onto white dwarfs we take into account the composition of 
accreted material (see §5.7). 

5.7. Mass Accumulation onto White Dwarfs 

A number of important phenomena, e.g., novae and 
Type la SN explosions or accretion-induced collapses, are 
associated with mass accretion onto WDs. We incorpo- 
rate the most recent results to estimate the accumulation 
efRciencies on WDs. In particular we consider accretion of 
matter of various compositions onto different WD types. 
We also include the possibility that NS formation can oc- 
cur via accretion induced collapse (AIC) of a massive ONe 
white dwarf (e.g., Bailyn & Grindlay 1990; Belczynski & 
Taam 2004a). 

In this section we discuss the accumulation of material 
and growth of the WD mass in binary systems. Only dur- 
ing dynamically stable RLOF phases can the mass accre- 
tion onto WDs be sustained for a prolonged period of time 
and hence affect the evolution of accreting WDs. During 
dynamically unstable cases (i.e., CE evolution) we assume 
that the WDs do not accrete any material. 

If dynamical instability is encountered for a binary with 
two white dwarfs we assume that a merger occurs. Based 
on the results of Saio & Nomoto (1998) mergers of a ONe 
WD with any type of WD companion and two CO WDs 
lead to either AIC and NS formation (if total merger mass 
A^morgcr IS abovc Afocs = 1.38 Mq) Or the formation of the 
single ONe WD (with new mass equal to Mmerger)- For 
mergers of CO WD and He WD, we assume a Type la SNa 
explosion; either sub-Chandrasekhar (Mmorgor < 1.44 Mq) 
(see Woosley, Taam, & Weaver 1986; Woosley & Weaver 
1994) or Chandrasekhar mass SN la (Mmergcr > 1-44 Mq. 
Mergers of other types of WDs have total mass below the 
Chandrasekhar mass, and in particular for CO WD and 
H WD wc assume formation of single CO WD, while for 
He WD and H WD we assume formation of single He WD 
with masses equal to Mmorgcr- 

During a phase of sustained mass accumulation the mas- 
sive ONe WD {K = 12) may eventually collapse to a 
NS. We include AIC in our standard model calculations 
since it naturally follows from the adopted accumulation 
physics (see below). Since little is known about potential 
asymmetries of the collapse, we either apply no natal kick 
(standard model) or a full natal kick (parameter studies) 
obtained from Arzoumanian, Chernoff, & Cordes (2002) 



or Hobbs ct al. (2005, sec also §6.2). However, we also 
allow for the possibility of SN la explosion instead of AIC 
in parameter studies. It is also worth noting the differ- 
ence between accretion and accumulation. The calcula- 
tion of accretion rate during stable RLOF was described 
in §5.1, and this rate could be used to calculate, for ex- 
ample, the accretion luminosity of the system (mostly in 
the UV part of spectrum for WD accrctors). However, it 
is believed that in many cases (see below) not all of the 
accreted material remains on the surface of the accreting 
WD. Mass is lost either in shell flashes (nova- like explo- 
sions) or through optically thick winds from the surface of 
accreting WDs. To calculate the actual WD mass growth 
through the RLOF phase the accumulation efhciency, Tyacuj 
which is defined as 

Macu = ?7acuA4cc (59) 

must be known. Here, Macu is the mass accumulation rate 
on the surface of WD and the mass accretion rate (Mace) 
is given by eg . ( [50)1 . In what follows we discuss the accu- 
mulation efficiency in various evolutionary scenarios. 

Accretion onto Helium and Hybrid white dwarfs. It is 
assumed that if the mass accretion rate Mace from the 
H-rich donor (ifdon = 0,1,2,3,4,5,6,16) is smaller than 
some critical value Mcriti, an unstable hydrogen shell flash 
will occur in the accreted layer. In response, the envelope 
will expand beyond the Roche lobe of the white dwarf. We 
shall assume no material is accumulated, and that the ac- 
cumulation efficiency is 77acu = 0.0, i.e. the entire accreted 
material is lost from the binary. If Mace > Mditi then the 
material piles up on the WD leading to mass loss from the 
system. Assuming that a contact binary configuration is 
not formed in this case, the system will eventually undergo 
an inspiral. For giant-like donors we assume the system 
evolves through CE to examine if the system survives; for 
all other donors we call it a merger and halt binary evolu- 
tion. The critical accretion rate is calculated as 

A/criti - loM^cciX * Q)-' Mq jt-' (60) 

where, Q = 6 x 10^^ erg g~^ is an energy yield of hy- 
drogen burning, X is the hydrogen content of accreted 
material. For Population I stars (metallicity Z > 0.01) we 
use X 0.7, ^0 = 1995262.3, A = 8, while for Population 
II stars {Z < 0.01) we use X = O.SJo = 31622.8, A = 5 
(Ritter 1999, see his eq. 10,12 and Table 2). 

If the mass accretion rate from the He-rich donor 
{Kdon = 7,8,9,10,17) is higher than Afcrit2 = 2 x 
10^* M0 yr~^ all the material is accumulated (?7acu = 1-0) 
until the accreted layer of material ignites in a helium shell 
flash. At this point degeneracy is lifted, a main sequence 
helium star (i^ace = 7) is formed and further accretion 
on the helium star is then taken into account. Following 
the calculations of Saio & Nomoto (1998) we estimate the 
maximum mass of the accreted shell at which the flash 
occurs as 

^j^j- ^ f -7.8 X lO^Afaee + 0.13 Afacc < 1-64 X 10"^ 

\ (instantaneous flash) Aface > 1.64 x 10~^ 

(61) 

where Mace is expressed in Mq yr ^. 

The newly formed helium star may overfill its Roche 
lobe, in which case either a single helium star is formed (He 
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or Hyb WD companion, A'don = 10, 17), a helium contact 
binary is formed (HeMS companion, A'don = 7) which we 
assume leads to a merger or the system goes through CE 
evolution (evolved helium star companion, if don = 8, jfl)- 
For accretion rates lower than Mcrit2, accumulation is 
also fully efficient (?7acu = 1-0). However, the SN la occurs 
at a sub-Chandrasekhar mass 

AfsNIa = -4 X lO^Macc + 1.34 Mq , (62) 

where Mace is expressed in Mq yr~^. For mass accretion 
rates close to Mcrit2, the above extrapolations from the 
results of Hashimoto et al. (1986) yield masses smaller 
than the current mass of the accretor, and we assume an 
instantaneous SN la explosion. We note that above ex- 
plosions disrupt the accreting WD, and although possibly 
subluminous, they appear as Type la SNe (no Hydrogen). 
We do not consider the accumulation of heavier elements 
since they could only originate from more massive WDs 
(e.g., CO or ONe WDs), which would have smaller radii 
and could not be donors to lighter He or Hyb WDs. 

Accretion onto Carbon/Oxygen white dwarfs. In this 
case we adopt the prescription from Ivanova & Taam 
(2004). For H-rich donors and mass accretion rates lower 
than 10~^^ Mq yr~^ there are strong nova explosions and 
no material is accumulated (?7acu = 0.0). In the range 
10~^^ < Mace < 10~® Mq yr~^ we interpolate for 77acu 
from Prialnik & Kovetz (1995, see their Table 1). For rates 
higher than 10"^ Mq yr~^ all accreted material burns into 
helium (?/acu = 1.0). Additionally we account for the 
effects of strong optically thick winds (Hachisu, Kato & 
Nomoto 1999), which eject any material accreted over the 
critical rate 

Merita = 0.75 X 10-'^(Macc - 0.4) Mq yr-\ (63) 

This corresponds to Jjacu = Mcrits/Macc for Mace > Merita- 
The accretor is allowed to increase in mass up to 1.4 Mq, 
and then explodes as a Chandrasekhar mass SN la. In the 
case of He-rich donors, if the mass accretion rate is higher 
than Aferit4 helium burning is stable and contributes to 
the accretor mass (77acu = 1.0). For rates in the range 
Meiit4 Merits accumulation is calculated from 

rfj, = -0.35(log Afaee + 6.1)^ + 1.02 [-6.5 ~ -6.34] 

^acu = -0.35(logMaec + 5.6)^ + 1.07 [-6.88 ^ -6.05] 

vlcn = -0.35(logA/aec + 5.6)^ + 1.01 [-6.92^-5.93] 

1.1-1.2^/ 0.541ogMaee + 4.16 [-7.06^-5.95] 

\ -0.54(logMaee + 5.6)2 + 1.01 [-5.95 -5.76] 

'?aeu = -0.175(logA/acc + 5.35)^ + 1.03 [-7.35^-5.83] 

rjl;f = -0.115(logMaee+5.7)2 + 1.01 [-7.4+-6.05] (64) 

and represents the amount of material that is left on the 
surface of the accreting WD of a specific mass (denoted 
by a superscript on Tjacu in Mq) after the helium shell 
flash cycle (Kato & Hachisu 1999, 2004). Logarithms of 
critical mass accretion rates for a given specific WD mass 

^For Xdon = 8 merger is assumed in such a case. 



are given in square brackets: [log(Mcrit5/ Mq yr ^) -f- 
log(Mcrit4/ Mq yr~^)]. To obtain the accumulation rate 
for CO WD within the mass range 0.7 — 1.4 Mq we in- 
corporate results of the closest (by mass) model from the 
set of cqs. [Ml If the WD mass drops below 0.7 Mq we 

use 77acu = 1-0 and we set logMerit4 = log Merits = -7.6 
(see Kato & Hachisu 2004). The mass of the CO WD 
accretor is allowed to increase up to 1.4 Mq, and then a 
Chandrasekhar mass SN la takes place in the two above 
He-rich accretion regimes. If mass accretion rates drop 
below Merits, the helium accumulates (?7acu = 1-0) on top 
of the CO WD and once the accumulated mass reaches 
0.1 Mq (Kato & Hachisu 1999), a detonation follows and 
ignites the CO core leading to the disruption of the ac- 
cretor in a sub-Chandrasekhar mass SN la (e.g., Taam 
1980; Garcia-Senz, Bravo & Woosley 1999). If the mass of 
the accreting WD has reached 1.4 Mq before the accretion 
layer has reached 0.1 Mq then the accretor explodes in a 
Chandrasekhar mass SN la. Carbon/Oxygen accumula- 
tion takes place without mass loss (r^acu = 1-0) and leads 
to SN la if Chandrasekhar mass is reached. 

Accretion onto Oxygen/Neon/Magnesium white dwarf. 
Accumulation onto ONe WDs is treated the same way as 
for CO WD accretors. The only difference arises when an 
accretor reaches the Chandrasekhar mass. In the case of 
ONe WD this leads to an AIC and NS formation, and bi- 
nary evolution continues (see Belczynski & Taam 2004a, 
2004b). 

6. SPATIAL VELOCITIES 

6.1. Overview 

All stars (single and binary systems) may be initialized 
with arbitrary velocities appropriate for a given environ- 
ment. For example, a galactic rotation curve may be used 
for a field population of a given galaxy, or a velocity disper- 
sion can be applied for a cluster population. The velocities 
of stars are then followed throughout their evolution. Sin- 
gle stars and binary systems are subject to recoil (change 
of spatial velocity) in SN explosions. Additionally, binary 
systems may be disrupted as a result of an especially vio- 
lent explosion. We account for both mass/angular momen- 
tum losses as well as for SN asymmetries (through natal 
kicks that NSs and BHs receive at their formation; see be- 
low). The detailed description of SN explosion treatment 
is given in BKB02. Here, we only list the new additions to 
StarTrack. The most important modification allows us to 
trace velocities of disrupted binary components after a SN 
explosion. For the first time, a full general approach with 
explosions taking place on orbits of arbitrary eccentric- 
ity (in contrast to circular orbits only) is applied to follow 
the trajectories of disrupted components. First population 
synthesis results are presented in Belczynski et al. (2006). 

6.2. Natal Kick Distribution 

At the time of birth, NSs and possibly BHs receive a 
natal kick, which is connected to asymmetries in SN ex- 
plosions. We use the distributions inferred from observed 
velocities of radio pulsars. We have replaced the natal 
kick distribution used in BKB02 (Cordes & Chernoff 1998) 
with two more recent alternatives. One presented by Ar- 
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zounmnian ct al. (2002) is a bimodal distribution with 
a weighted sum of two Gaussians, one with a ~ 90 km 
sec-i (40%) and another with a = 500 km sec^^ (60%). 
The other was derived by Hobbs et al. (2005) and is a 
single Maxwellian with a ~ 265 km sec~^. According to 
this most recent study there is no indication of a bimodal 
(low- and high-velocity) kick distribution claimed in ear- 
lier studies (e.g.. Fryer, Burrows & Benz 1998; Cordes & 
Chernoff 1998; Arzoumanian et al. 2002). If this is in- 
deed true, then some theoretical models built in support 
of the bimodal kick distribution (e.g., Pfahl et al. 2002b 
and Podsiadlowski et al. 2004 model of high mass X-ray 
binaries) may need to be revised. Motions of many hun- 
dreds of pulsars are expected to be measured in the next 
few years. These measurements will provide better con- 
straints on the natal kick distribution (Hobbs et al. 2005). 
Until then we will use both distributions to assess the as- 
sociated uncertainties in StarTrack calculations. 

Compact objects formed without any fall back receive 
full kicks drawn from one of the above distributions; this 
case includes most of NSs (sec §2.3.1). The only excep- 
tion arc NS formed through ECS, for which we adopt ei- 
ther no natal kicks (standard model) or full kicks (as part 
of parameter studies). During accretion induced collapse 
of WD to NS in an accreting binary system, the NS is 
formed through the same ECS process and the same pre- 
scription is adopted. The recent numerical simulations of 
AIC and NS formation through ECS point towards sig- 
nificantly lower energies of explosion than for regular core 
collapse SNe (Dessart et al. 2006), and although the kick 
mechanism is not yet identified, these results may be also 
indicative of lower kicks in ECS NS formation. For com- 
pact objects formed with partial fallback, heavy NSs and 
light BHs, kicks are lowered proportionally to the amount 
of fallback associated with NS/BH formation 



Mcick = (1 - ffh)V 



(65) 



where V is the kick magnitude drawn from cither Arzou- 
manian ct al. (2002) or Hobbs et al. (2005) distribution, 
and /fb is a fallback parameter, i.e., the fraction (from 
to 1) of the stellar envelope that falls back (sec also 
§2.3.1). For the most massive BHs, formed silently (no 
SN explosion) in a direct collapse (/fb = 1) of a massive 
star to a BH, we assume that no natal kick is imparted. 
The adopted natal kick distribution and kick scaling for 
NSs and BHs can be readily changed for parameter stud- 
ies (e.g., full BH kicks). 

6.3. Supernova disruptions 

Just prior to the SN explosion, the two components of 
the binary move with velocities v{ and tT|, which, in the 
center of mass (CM) system of coordinates, denoted here 
with the superscript /, satisfy 



(66) 



where Mi denotes the SN component and M2 its compan- 
ion. Subscripts int, fin stand for initial and final values. 

Wc make no assumptions about the orbit; it can have 
an arbitrary eccentricity, in contrast to the derivation by 
Tauris & Takcns (1998), who assumed that the orbit is cir- 
cular prior to the explosion. At the moment of a supernova 



explosion the orbital separation is rQU. The exploding star 
loses its envelope, its mass becomes Mi fin and receives a 
kick w, so now its velocity in the coordinate system / is 



-'l.int 



(67) 



The secondary star may be affected by the expanding 
shell and may receive an additional velocity Himp, how- 
ever it has been shown (Kalogera 1996) that the effect 
of this velocity is small, unless the pre-supernova orbital 
separation is smaller than ~ 3R0. Wc assume that the 
velocity of the companion is not affected by the impact of 
supernova ejecta. We also assume that the velocity of the 
shell is large and the shell leaves the system quickly, i.e. 
Wshoii >> ToP, where P is the orbital period of the system 
prior to the explosion. 

In order to calculate the final velocity of the two stars 
we first transform the velocities to the CM system of the 
two post-SN stars. The velocity of this system, denoted 
as //, in relation to system / is 



Mi,fi„w( int + M2vi 



Ml fin + M2 



(68) 



The relative velocity of the two stars in this system is 



^11 



Vn + W — Vu 



(69) 



while the initial direction between the two stars remains 



the same as in the coordinate system /, n 
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In this 



new system the relative motion of the stars is a hyper- 
bola in the plane perpendicular to the angular momentum 
vector: 

(70) 



J = /iro^^^ X v'^ 



where /i = Mi3nM2/(Mi^fin -I- M2) is the reduced mass 
of the system. It is convenient now to introduce a third 
coordinate system /// in which the angular momentum J 
lies along the z-axis. The transformation from // to III 
is a rotation TZ: v^^^ = TZv^' ,n^^^ = TZn^' . The orbit in 
/// is described by 



P 



1 + e cos ( 



where 



(71) 



(72) 



with E = iJi{v^^Y /2 — a/|ro| is the (positive) energy of 
the system, and a = GAjfi_flnM2. The final velocity, at 
r ^ 00, follows from energy conservation: 





(73) 



In order to find the direction of the final velocity we 
note that conservation of angular momentum implies that 
at infinity (r 00): the final relative is parallel to 
the direction between the stars J^jf,^. The initial position 
of the two stars on the trajectory described by eg. [TTlis 



1 I p 

cos .^int = - 



?'0 



(74) 
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The sign of the angle ipi^t is negative if the two stars 
initially lie on the descending branch of the hyperbola 



^///fg^^ > and positive if they are on the ascending 

one "uf/Z^o^^ < 0- ^^^'^ ^^^^ case, when the two stars are 
initially on the descending branch, we need to compare the 
distance of closest approach on the orbit rmin = + e) 
with the radius of the companion star to examine whether 
the two stars collide instead of escaping to infinity. 
We obtain the final position on the trajectory from 



1 



cos (/Jfin 



(75) 



and (/3fin > 0. Thus the final direction between the two 



stars at r = oo is flf^l^^ = T{(pfin — fint)n'^^ , where T((/)) is 
the matrix of rotation around the z-axis, and their relative 
velocity is: 




III 



(76) 



We now have to transform ciuantities from system /// back 
to system / to obtain the final velocities of the two dis- 
rupted binary components in the initial (pre-SN) CM sys- 
tem: 



= n 



-1 / 



"cm- 



(77) 
(78) 



7. DISTRIBUTIONS OF INITIAL PARAMETERS 



Each binary system is initialized by four parameters, 
which are assumed to be independent: the primary mass 
Ml (the initially more massive component), the mass ratio 
q = M2/M1, where M2 is the mass of the secondary, the 
semi-major axis a of the orbit, and the orbital eccentricity 
e. 

For both single stars and binary system primaries, we 
use the initial mass function adopted from Kroupa, Tout 
& Gilmore (1993) and Kroupa & Weidner (2003), 



«'(Mi) cx 




0.08 < Ml < 0.5 Mq 
0.5 < Ml < 1.0 Mq 
1.0 < Ml < 150 Mq 



(79) 



where parameter Ofimf = 2.35 — 3.2, with our standard 
choice being 2.7 for field populations and 2.35 for cluster 
populations. Stars are generated within an initial mass 
range: Mmin — Mmax, and the range is chosen accordingly 
based on the targeted stellar population. For example, NS 
studies would require evolution of single stars within range 
~ 8 — 25 Mq while the formation of WDs would require an 
initial range ^ 0.08 — 8 Mq. Binary evolution, due to mass 
transfer events (both mass accretion and mass loss) may 
significantly broaden any of the ranges mentioned above. 
We assume a flat mass ratio distribution. 



(80) 



in the range g = — 1 in agreement with recent obser- 
vational results of Kobulnicki, Fryer & Kiminki (2006). 
However, it should be noted that massive binaries may 



form with components of comparable mass (Pinsonneault 
& Stanek 2006), and we will test this alternative mass ra- 
tio distribution in our parameter studies. Given the value 
of the primary mass and the mass ratio, we obtain the 
mass of the secondary M2 = qM 1 . 

The distribution of initial binary separations is assumed 
to be flat in the logarithm (Abt 1983), 



r(a) cx — , 
a 



(81) 



where a ranges from a minimum value, such that the pri- 
mary fills at most 50% of its Roche lobe at ZAMS, up to 
10^ Rq. 

Finally, we adopt the thermal-equilibrium eccentricity 
distribution for initial binaries, 



S(e) = 2e, 



(82) 



in the range e 
Mayor 1991). 



0—1 (e.g., Hcggic 1975; Duqucnnoy & 



8. CALIBRATIONS AND COMPARISONS 

8.1. Mass Transfer Sequences 

In the following subsections we present Star Track mass 
transfer calculations and compare them to published and 
unpublished results based on the use of stellar evolution 
and mass transfer codes. 

8.1.1. Case B Mass Transfer: MS+HG binary 

We choose this RLOF sequence from Wellstein et al. 
(2001, their Model B) and start with a 16 Mq -f 15 Mq 
ZAMS binary in a 8 day circular orbit. RLOF starts af- 
ter the primary evolves off the MS. The system at the 
onset of RLOF (t = 11.5 Myr since ZAMS) is character- 
ized by: Ki = 2, K2 = 1, Porb = 7'^.9, e = 0,Mi = 
15.6 Mq, M2 = 14.7 Mq, Ri = 19.9 Rq, and R2 = 
11 A Rq. The evolution of the system during the RLOF 
phase is shown in Figure [S] 

The RLOF phase proceeds on the thermal timescale of 
the donor, which is rapidly expanding while crossing the 
Hertzsprung gap. First, there is a phase characterized by 
very high mass transfer rates (~ 3 x 10"'^ Mq yr~^), until 
the mass ratio is reversed and the donor becomes the less 
massive binary component. Shortly thereafter, the trans- 
fer rate slowly decreases (~ 10^"^ — 10^'' Mq yr^^). The 
mass ratio is reversed just right after the onset of RLOF 
and the orbit expands. 

Our assumption is that all (100%) of the transferred 
material is accreted by the companion (conservative evo- 
lution is adopted for a better comparison with Wellstein 
et al. 2001 calculation). RLOF terminates when the en- 
velope of the donor is nearly exhausted and its radius 
contracts below the Roche lobe radius, thereby, caus- 
ing the system to become detached. The primary loses 
most of its mass and becomes a core helium burning star 
(Ml = 3.99 Mq, Ki = 4), while the secondary gains mass 
and is rejuvenated (M2 = 26.2 Mq, K2 = 1). The orbital 
period increases to reach ~ 80 days at the end of the RLOF 
phase. Both stars continue to evolve in a detached config- 
uration. Eventually, the primary becomes a naked helium 
star {Ml = 3.98 Mq, Ki = 7) that evolves and loses some 
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mass {Ml — 3.5 Mq, Ki = 8 — 9) until the final explo- 
sion and forms a neutron star [Mi = 1.45 Mq, = 13). 
Since the mass of the primary is quite large during the he- 
lium burning stage it does not significantly expand and in 
particular it does not initiate another RLOF. The helium 
star primary reaches a maximum radius of i?i ~ 10 
just prior the explosion, while its Roche lobe radius is 
^i.iob ~ 60 Rq. After the explosion the system is dis- 
rupted due to a large natal kick the neutron star receives 
(e.g., Hobbs et al. 2005). At this point the secondary 
is still on the main sequence (M2 = 25.7 M©, K2 = 1), 
and will eventually (in ^ 2 Myr) form a single black hole 
(M2 = 6.4 Mq, Ki = 14). 

The calculation of Wellstein et al. (2001) shows simi- 
lar behavior during the first RLOF phase in terms of the 
duration, mass transfer rate, and orbital period. The fi- 
nal donor and accretor masses in both simulations are al- 
most the same. Therefore, our calibration and rates used 
for thermal and then nuclear timescale RLOF are in very 
good agreement with the detailed evolutionary calculation 
of Wellstein et al. (2001). 

However, Wellstein et al. (2001) find a second RLOF 
phase, when the primary expands again, once it becomes 
a helium giant {Ki = 8—9) after it loses significant amount 
of mass in a stellar wind: AIi = 3.8, 2.8 Mq at helium star 
formation and at the start of the second (case BB) mass 
transfer phase, respectively. Their system evolves through 
the, so called, case BB mass transfer phase. The RLOF 
stops when the primary loses a significant part of its he- 
lium envelope. The system becomes wider, and eventually 
the primary explodes in type Ib/c supernova forming a 
neutron star. The difference in the evolution into a second 
RLOF phase is explained by the difference in modeling the 
helium star evolution, and in particular the stellar winds. 
The helium star primaries in both calculations start with 
about the same mass (Mi ~ 4 Mq) and their lifetimes are 
similar (~ 1 Myr). However, in the case of Wellstein et al. 
(2001) the helium star loses significant amount of mass: 
~ 1 M0, while our star loses only half of that: ~ 0.5 Mq. 
Therefore, our helium star remains massive and does not 
significantly expand to initiate the second RLOF as it is 
found by Wellstein et al. (2001). The reason for this dif- 
ference in wind mass loss from helium stars is that we 
use the downward revised empirical wind mass loss rates 
(Hamann & Koesterke 1998; see also Hurley et al. 2000) 
that take into account wind clumping and predict rates 
lower by factor of ~ 2 than previously estimated (Hamann 
Schonberner, & Heber 19820. 

8.1.2. Case A Mass Transfer: MS+MS binary 

This RLOF sequence is selected from Wellstein et al. 
(2001, their Model A). We start with a 12 Mq 7.5 Mq 
ZAMS binary in a 2.5 day circular orbit. RLOF starts 
while the primary still evolves through the MS phase. The 
system at the onset of RLOF {t = 14.8 Myr since ZAMS) 
is characterized by: Ki ~ 1, A'2 = 1, Porb = 2''. 3, e = 
0, Ml = 11.9 Mq, M2 = 7.5 Mq, i?i = 8.3 R©, and R2 = 
4.0 Rq. The evolution of the system during the RLOF 
phase is shown in Figure [Sj In this calculation we invoke 
conservative evolution (/a = 1.0; all mass lost from donor 

^Wellstein et al. (2001) refer to Wellstein & Langer (1999) for the employed mass loss rates. In Wellstein & Langer (1999) we find that both 
old and revised mass loss rates for helium stars are presented. However, it is apparent from the difference in the results that Wellstein et al. 
(2001) evolution of their Model B binary employs the old (high) rates. 



is accreted by the companion) to match the assumption in 
Wellstein et al. (2001). 

First phase. At first, the RLOF proceeds on the thermal 
timescale with a mass transfer rate of ^ 5 x 10~^ Mq yr~^, 
through the so called rapid case A transfer phase. The 
transfer rate then rapidly decreases by more than 2 or- 
ders of magnitude until the component masses are nearly 
equal. Subsequent evolution proceeds on the much slower 
nuclear timescale of the donor with transfer rates below 
10~^ Mq yr~^. RLOF continues until the final stages of 
the donor MS lifetime, when the primary contracts and 
detaches from its Roche lobe. The evolution of the or- 
bital period is characterized by an initial small decrease 
and then (after the thermal timescale phase has ended) a 
slow but rather constant increase up to 3.5 days. At that 
point the primary mass is ~ 6 M© and the secondary mass 

- 13 Mq. 

Second phase. After ~ 0.5 Myr the primary starts 
expanding as it enters the Hertzsprung gap and RLOF 
restarts. This mass transfer phase is much more rapid 
and is driven by expansion of the primary. This phase 
is characterized by high mass transfer rates (10~^ — 
10~^ Mq yr~^) and the envelope of the primary is soon 
(~ 0.3 Myr) exhausted, ending the second RLOF phase. 
During this relatively short phase, the orbit expands signif- 
icantly (final orbital period ~ 260 days) , while the primary 
loses most of its mass {Mi sa 1 Mq) while the secondary, 
still on its MS, gains mass (M2 « 18 Mq) and is rejuve- 
nated. The dramatic orbit expansion is an effect of the 
rather extreme mass ratio for this system at the time of 
the second RLOF. For both RLOF phases conservative 
evolution was applied. The evolution of this system ends 
when the massive secondary evolves off MS, initiating a 
CE phase while crossing the Hertzsprung gap. This phase 
leads to inspiral and merger. 

The calculation of Wellstein et al. (2001) shows a 
qualitatively similar system behavior during both RLOF 
phases; initial high mass transfer phase, then a slower 
one, short break in RLOF followed by another rapid phase 
while the donor evolves off the MS. Also the mass transfer 
rates are comparable as it is duration of the second RLOF 
phase. 

However, there is a difference in the duration of the first 
RLOF phase, our calculation showing factor of 3 longer 
RLOF phases than that of Wellstein et al. (2001). The 
initial rapid (on a thermal timescale) phase of RLOF lasts 
longer in the Wellstein et al. (2001) calculation, and the 
mass transfer rates are slightly different than in our ap- 
proximations. These result in a somewhat different mass 
evolution of the binary components (e.g., our secondary 
has reached 13.5 M© at the end of first RLOF phase while 
in Wellstein et al. (2001) calculation it ends up with 
5.5 Mq), and therefore different mass ratio of the system 
that in turn influences the subsequent mass transfer cal- 
culations and period evolution. Additionally, we include 
spin-orbit coupling in our calculations. The period of our 
system after the second RLOF is longer (260 days) than 

- 100 days found by Wellstein et al. (2001). 

We note that Wellstein et al. (2001) account for the 
rejuvenation effects in detail, given that they use a stel- 
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lar structure code and do not have to assume full rejuve- 
nation as it is adopted by Tout et al. (1997), Hurley at 
al. (2002), and our code implementation. Wellstcin et al. 
(2001) remark that had full rejuvenation been assumed in 
their calculation the system would have ended in a CE 
merger during the expansion phase of the secondary after 
MS evolution in agreement with our findings. 

8.1.3. BH-MS binary 

This calculation starts with a 10 M© BH + 5 Mq ZAMS 
star. We let the secondary evolve through about half of 
its MS lifetime before bringing the system into contact at 
t = 51.3 Myr (counted from the secondary ZAMS) . The 
system at the onset of RLOF is characterized by: Ki = 
14, K2 = I, Porb = I'^.O, e = 0, Ml = 10 Mq, M2 = 
5Mq, Ri = 0.000042 Rg, and i?2 = 3.5 Rq. The evo- 
lution of the system during the RLOF phase is shown in 
Figure El 

First phase. RLOF is stable and proceeds on the nu- 
clear timescale of the secondary with a mass transfer rate 
of ~ 2 X 10~® Mq yr~^. Since this rate is sub-Eddington 
we allow all the transferred material to be accreted onto 
the primary BH, which increases its mass to ~ 11.5 Mq, 
while the secondary mass decreases to ~ 3.5 Mq. During 
this phase, the period increases from 1 to 2 days. The 
phase ends when the secondary begins contraction at the 
end of its MS life. 

Second phase. RLOF restarts when the secondary 
crosses the Hertzsprung gap with mass transfer proceed- 
ing at the high rate (~ 10~^ Mq yr~^) corresponding to 
rapid expansion of the star on its thermal timescale dur- 
ing that phase. At some point the donor starts ascending 
along the red giant branch, and the transfer rate drops 
by about an order of magnitude to 3 x 10"'' Mq yr~^. 
Since the transfer rate is super-Eddington throughout this 
entire phase we limit accretion onto the BH to the Ed- 
dington rate, allowing the rest of the material to leave the 
binary with the specific orbital angular momentum of the 
BH. In the end the BH has increased its mass to 12.6 M© 
and the mass of the donor has decreased to 0.6 M©. The 
orbit expands significantly (^ 300 days) during this rapid 
RLOF phase. 

The RLOF phase ends at the point when the donor, due 
to the loss of its almost entire H-rich envelope, stops its 
expansion. The system ends its life as a wide BH-WD 
binary. 

The same RLOF sequence was calculated with the de- 
tailed stellar evolution code of Ivanova et al. (2003; also 
see Ivanova & Taam 2004). The comparison of the two 
phases of RLOF shows overall qualitative agreement with 
the StarTrack calculation. The mass transfer rates are 
virtually the same: ~ 2 x 10~* and ~ 10~^ M© yr~-^, for 
first and second phase, respectively. However, the detailed 
calculation with the evolution code shows a longer dura- 
tion (by a factor ^ 2) for the first RLOF phase. 

8.1.4. BH-RG binary 

This calculation starts with a 7 M© BH 2 M© ZAMS 
star. We let the secondary evolve through about one third 
of its red giant lifetime before bringing the system into 
contact at t = 1180.4 Myr (counted from the secondary 
ZAMS). The system at the onset of RLOF is character- 
ized by: Ki = 14, K2 = 3, Porb = 4''.8, e = 0, Mi = 



7M©, M2 = 2M©, Ri = 0.000030 R©, and R2 = 
7.1 R©. The evolution of the system during RLOF phase 
is shown in Figure [3 

RLOF is stable and proceeds through the entire RG 
phase {K2 = 3) on the nuclear timescale of the donor. The 
mass transfer rate is sub-Eddington and thus the material 
transferred to the BH is entirely accreted. In the end the 
mass of the BH is increased to 8.4 M© while the mass of 
the donor is decreased to 0.6 M©. As the donor expands, 
ascending the RG branch, the orbit expands as well, and 
finally the RLOF phase terminates at an orbital period of 
~ 90 days. The phase ends when the donor contracts upon 
igniting helium in its core. The system eventually forms a 
wide BH-WD binary. 

This RLOF sequence was also calculated with the de- 
tailed stellar evolution code of Ivanova et al. (2003). The 
mass transfer rates found in both cases are similar (~ 
10"''' - 10"® M© yr~i) and in this case the StarTrack 
timescales are shorter, but do not differ by more than 50%. 

8.1.5. Short period NS-RG binary 

This RLOF sequence is chosen from Tauris & Savonije 
(1999, their example 2b). We start with a 1.3 M© NS 
+ 1.6 M© ZAMS star in a 3 day circular orbit. RLOF 
starts while the secondary is on the RG branch (t = 2321.4 
Myr since secondary ZAMS) and the binary is described 
by: Ki = 13, K2 = 3, Porb = 2'*.8, e = 0, Mi = 
1.3 M©, M2 = 1.6 M©, Pi = 0.000014 R©, and R2 = 
4.7 R©. The evolution of the system during the RLOF 
phase is shown in Figure O 

At first the RLOF proceeds on a thermal timescale 
with a highly super-Eddington mass transfer rate (^ 
10"^ M© yr~^). After the donor becomes less massive 
than its accretor, the mass transfer is driven by the ex- 
pansion of the red giant donor (on its nuclear timescale) 
at a much smaller rate of ~ 10~® M© yr~^. As the mass 
transfer rate decreases, the NS starts to accrete efficiently 
and its mass increases to 1.9 M©. Eventually, after ~ 65 
Myr of RLOF, the RG secondary loses most of its mass 
(Af2 = 0.28 M©) and contracts, leaving a remnant helium 
WD. At this point the RLOF phase ends (orbital period 
60 days), and further evolution leads to the formation of 
wide binary, with a recycled pulsar. 

Comparison with the detailed evolutionary calculation 
of Tauris & Savonije (1999) shows good agreement be- 
tween the results. The detailed calculations show an ini- 
tial rapid RLOF phase followed by a sub-Eddington mass 
transfer phase, eventually leading to the formation of NS- 
He WD binary. Final component masses (NS and donor: 
2.05 and 0.29 M©, respectively) are very similar to the 
ones obtained with StarTrack. The final orbital period of 
42 days obtained by Tauris & Savonije (1999) is shorter 
than in our calculation (60 days). In addition, there is 
a difference in the duration of RLOF phase, lasting 123 
Myr in the Tauris & Savonije (1999) model, as compared 
to 60 Myr in our calculations. This may be understood in 
terms of a different treatment of binary interactions (tides, 
magnetic braking, winds) as well as the difference in stel- 
lar models which may lead to a different starting point of 
RLOF. 

8.1.6. Long period NS-RG binary 
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This RLOF sequence is taken from Tauris & Savonije 
(1999, their example 2c). We start with a 1.3 Mq NS + 
1.0 Mq ZAMS star in a 60 day circular orbit. RLOF starts 
while the secondary is on the RG branch (t = 12312.5 
Myr since secondary ZAMS) and the binary is described 
by: Ki = 13, K2 = 3, Porb = 60'^.708033, e = 0, Mi = 
1.3 Mq, M2 = 0.98 Mq, Ri = 0.000014 Rq, and R2 = 
30.5 Rq. The evolution of the system during RLOF phase 
is shown in Figure [TUl 

RLOF is highly super-Eddington and driven by the ex- 
pansion of the donor on a nuclear timescale. Only shortly 
before the system detaches as a result of the exhaustion 
of the donor's envelope, the transfer rate becomes sub- 
Eddington. As a result, the donor loses most of its mass 
(M2 = 0.4 M0) while the NS hardly accretes any ma- 
terial (Ml = 1.43 Mq). The orbit expands throughout 
this phase with the orbital period increasing to over 300 
days. The system eventually forms a wide NS-He WD bi- 
nary, with a potential recycled pulsar (the NS has accreted 
~ 0.1 Mq). 

The above results are very similar to the calculations of 
Tauris & Savonije (1999), who obtain a 1.5 Mq NS with a 
0.4 Mq NS-He WD binary in a 382 day orbit. The mass 
transfer rates and duration of the RLOF phases arc similar 
in both calculations. 

8.1.7. Long period NS-He star binary 

This RLOF sequence follows from Dcwi & Pols (2003, 
see their Fig. 1). We start with a 1.4 Mq NS + 2.8 M© 
ZAMS He star in a 10 day circular orbit. RLOF starts 
while the secondary is already an evolved He star {t = 2.9 
Myr) and the binary is described by: Ki = 13, K2 = 
8, Porb = 9'*.804617, e = 0, Mi = 1.4 M©, AI2 = 
2.5 Mq, Pi = 0.000014 Rq, and P2 = 15.4 Rq. The 
evolution of the system during RLOF phase is shown in 
Figure [TlJ 

RLOF proceeds on the donor's thermal timescale 
throughout the entire phase. The very high mass transfer 
rate (6 x 10"'^ Mq yr~^) makes this phase very short and 
RLOF stops after the envelope of He star is exhausted. 
Since the mass transfer rate is highly super-Eddington. 
the NS hardly accretes any material while the donor loses 
its entire He-rich envelope (M2 = 1.65 Mq). The orbital 
period at first decreases to reach a minimum at 8.9 days, 
and then increases to 9.5 days at the end of RLOF phase. 
After the phase of RLOF the secondary core explodes in 
SN Ic leading to double neutron star formation (provided 
that a natal kick does not disrupt the binary). This result 
was presented also in Ivanova et al. (2003). 

Dcwi & Pols (2003) calculated mass transfer rates span- 
ning the range: 10~^— 10~^ Mq yr~^. Our rate is constant 
and close to the high end of the Dewi & Pols (2003) range. 
We have adopted a constant mass transfer rate following 
Paczynski (1971) who pointed out that thermal timescale 
mass transfer rates do not vary by more than factor of 2-3 
(for details see § 5.3). This system may appear as an X-ray 
binary during this phase. However, the chances of catching 
it at this phase are very small, since the thermal timescale 
mass transfer is very short. Besides, in this case the X- 
rays may be significantly degraded because of high optical 
depths (material shed out of the system). On the other 
hand, some of these sources might appear to be soft 7-ray 
emitters (i.e. 20 — 100 keV range, tail of X-ray emission) 



with high intrinsic absorption, and the discovery of ob- 
jects with these broad characteristics (see e.g.. Dean et al. 
2005) lends some hope for detecting this phase of binary 
evolution. The results from Dewi & Pols (2003) reveal a 
different period evolution than in our simulation; RLOF 
starts at higher value (10.46 days), and then decreases to 
10.37 days. However, the period changes in both calcu- 
lations are rather small, and arc probably related to our 
consistently high mass transfer rate throughout the RLOF 
phase. This leads to higher mass and angular momentum 
loss from the binary which determines the orbit evolution. 
Additionally, we include tidal interactions between binary 
components (see §3.3 and §8.2 ). These differences be- 
tween models for low mass helium stars were already noted 
by Dewi & Pols (2003). 

8.1.8. Short period NS-He star binary 

We choose this RLOF sequence from Dewi & Pols 
(2003, see their Fig. 3). We start with a 1.4 Mq NS 
3.6 Mq ZAMS He star in a 0.6 day circular orbit. RLOF 
starts while the secondary is already an evolved He star 
{t = 2.0 Myr) and the binary is described by: Ki = 
13, A'2 = 8, Poib = 0^^.59, e = 0, Mi 1.4 Mq, M2 = 
3.2 Mq, Pi = 0.000014 Rq, and P2 = 2.4 Rq. The evo- 
lution of the system during RLOF phase is shown in Fig- 
ure [TH 

The RLOF proceeds on the donor's thermal timescale 
with a mass transfer rate of ~ 10"'^ Mq yr~^ until the 
envelope of He star is almost exhausted. Since the mass 
transfer rate is highly super-Eddington, the NS does not 
accrete much material while the donor loses most of its He- 
rich envelope (M2 = 2 Mq). The orbital period decreases 
from 0.6 to ~ 0.4 days at the end of RLOF phase. Af- 
ter the RLOF phase ceases the secondary explodes in SN 
Ib/Ic leading to a close double neutron star system (again 
provided that a natal kick does not disrupt the binary). 

The Dewi & Pols (2003) RLOF sequence for this case 
is very similar to our calculation. They find a period de- 
crease (from 0.65 to 0.47 days) and a high constant mass 
transfer rate of a few xlO""* Mq yr~^. The inspiral phase 
and CE is not expected in this case, and therefore further 
evolution may lead to a close double neutron star forma- 
tion. 

8.2. Tidal Evolution Calibration 

Whenever coeval binary populations in nearby clus- 
ters are observed to constrain the circularization rate, it 
is found that standard tidal dissipation theories do not 
match the data (see Meibom & Mathieu 2005 for a recent 
review). In all cases an increase in the tidal dissipation 
rate appears necessary (Claret & Cunha 1997; Terquem et 
al. 1998). Depending on which theory is used, the increase 
needed in the overall efficiency of tidal dissipation is by a 
factor ~ 10 - 100. 

We have used StarTrack models to calibrate our theo- 
retical treatment by comparing them against observations 
of (i) the cutoff period for circularization in a population 
of MS binaries (in M67), and (ii) the orbital decay accom- 
panying tidal synchronization in a high mass X-ray binary 
(LMC X-4). The results, presented in two following sub- 
sections, confirm that tidal dissipation, at least in case of 
convective stars, is more effective than predicted by our 
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simple theory. Therefore, in all our standard model calcu- 
lations, we will use an increased rate of tidal dissipation 
for convective stars, corresponding to -Ftid,con = 50, while 
using standard dissipation for radiative stars -Ftid,rad = I7 
but we will also allow for even more effective tidal dis- 
sipation rates in our parameter studies (all the way to 
^tid,con = 100 and i^tid.rad = 100). See §3.3 for our imple- 
mentation of tidal dissipation theory and the definition of 

-Ptid- 

8.2.1. Cutoff Period for M67 

Open star clusters have often been used to test tidal in- 
teraction theories (Mathieu et al. 1992; Meibom & Math- 
ieu 2005). Observations of single- and double- line spectro- 
scopic binaries allow estimates of the periods and eccen- 
tricities for a number of systems within clusters. It was 
expected and then confirmed that the cutoff period (Pcut , 
the longest period of a circular binary) should increase 
with the age of the cluster. The tidal dissipation depends 
strongly on the orbital separation and therefore the wider, 
longer-period binaries will take a greater time to circular- 
ize. In principle, with knowledge of the initial conditions in 
a given cluster, the observed value of the cutoff period may 
be used to calibrate the efficiency of tidal interactions. In 
practice, binaries within a given cluster form with eccen- 
tricities, separations and angular momenta which are not 
precisely known. In addition, the observed samples may 
suffer from small number statistics (the observed cutoff pe- 
riods are only lower limits), rendering such a calibration 
quite uncertain. However, we can use the cutoff-period 
observations to provide at least an order of magnitude es- 
timate for the factor by which any standard theoretical 
estimate must be increased. 

M67 is an old open cluster with an age of 3.98 Gyr and 
observed cutoff period of 10 — 12 d (Mathieu, Latham & 
Griffin 1990; Mathieu et al. 1992) and a solar metallic- 
ity stellar population (Janes & Phels 1994). The period 
was estimated for a sample of MS binaries with compo- 
nents close to the cluster turnoff mass. Recently Meibom 
& Mathieu (2005) proposed a new way to estimate the 
point of transition from circular to eccentric systems. In- 
stead of a simple cutoff period, they use a new estima- 
tor called the "tidal circularization period." This period 
is found from fitting a special function which mimics the 
tidal circularization isochrone of the most frequently oc- 
curring eccentric binary orbits for a given cluster. They 
find that the tidal circularization period for M67 is 12.1 d. 

Several calculations, with different efficiencies of tidal 
dissipation, were performed to try to reproduce the bi- 
nary population of the open cluster M67. In each cal- 
culation we have evolved lO'' binaries at solar metallic- 
ity with component masses in the 0.7 — 1.4 Mq range, 
requiring that the mass ratio be greater than 0.5. The 
limits are somewhat arbitrary, but chosen to include the 
population of bright MS stars observed in M67. Most of 
these stars have convective envelopes, and therefore we 
try to calibrate the scaling factor for convective envelopes 
(-Ftid.con) while keeping the one for radiative stars constant 
(_F'tid,rad = 1)- Thc initial distributions were chosen as in 
our standard evolutionary model (see § 5.7), but with IMF 
exponent aimf = 2.35, which is more appropriate for clus- 
ters (Kroupa & Weidner 2003). 

In Figure [T3] we show synthetic binary MS popula- 



tions in the period-eccentricity plane corresponding to an 
evolution with different efficiencies for the tidal interac- 
tion. As expected we see that the cutoff period increases 
for more efficient tidal interactions, Pcut — 4, 7,10 d for 
^tid.con = 1, 10, 100, respectively. It is found that only for 
significantly increased dissipation (Ftid,con ^ 10 — 100) the 
the predicted cutoff period approach the observed value of 
10-12 days. An additional calculation with Ftid,con = 1000 
results in a cutoff period of ~ 16 days, now clearly higher 
than the observed value. 

8.2.2. Orbital decay of LMC X-4 

Levine, Rappaport & Zojcheski (2000) measured an or- 
bital period decay for the high mass X-ray binary (HMXB) 
LMC X-4. The system consists of a 1.3 Mg NS and a 
massive 15.6 Mq companion in a 1.4-day circular or al- 
most circular orbit (Woo et al. 1996; van der Meer et al. 
2005). The X-ray emission in HMXBs is believed to arise 
from wind accretion onto the compact object; however it 
was also suggested that some systems may be in an at- 
mospheric RLOF phase (e.g., Kaper 2001). For wind- fed 
detached systems, the orbital decay may be directly con- 
nected to the tidal interaction of the HMXB components. 
The secondary is a massive star, and the source of tidal 
dissipation is radiative damping. Therefore, we use LMC 
X-4 to check the efficiency of tidal interactions for radiative 
envelopes (i^tid,rad)- The rotation of the massive compo- 
nent decreases with time as it expands during its evolution. 
On the other hand, the tidal forces act to synchronize the 
massive component, resulting in loss of orbital energy and 
angular momentum, i.e., decay of the orbit. 

If, in fact, LMC X-4 is a wind-fed system and not in 
RLOF, then the massive star must be smaller than its 
Roche lobe i?roche = 8 R©. A 15.6 Mq star exceeds that 
size, while still on MS, after about 10.5 Myr of evolution 
(from the ZAMS). Subsequent RLOF is dynamically un- 
stable (extreme mass ratio) and leads to a rapid merger 
of the binary components, terminating the HMXB phase. 
We perform a set of calculations for a synthetic binary sim- 
ilar to the LMC X-4 using our standard model parameters, 
with a metallicity appropriate for the LMC {Z = 0.007). 
We assume that the binary configuration is detached and 
we calculate the rate of orbital decay. The orbital decay 
rate depends crucially on the current relative radius of the 
massive component of LMC X-4 (ex (R/a)^, see eg. [T5)) . 

The radius of the 15.6 Mq star {Z = 0.007) increases 
from i?2 = 4.5 Rq on ZAMS to R2 = 13 Rq at the end 
of MS phase, which takes ~ 13 Myr. Based on the Roche 
lobe radius of secondary for 1.4 day orbit the secondary 
fills its Roche lobe at 11 Myr, which is close to the 
end of MS phase. Since the primary has already evolved 
and has formed a NS, a significant amount of time must 
have elapsed since the binary formation. For example a 
30 — 35 Mq star takes ~ 6 Myr to form compact ob- 
ject, and such a massive star would have formed a NS 
only if stripped of a significant part of its mass in RLOF 
episode. For more massive primaries, the evolution would 
be slightly faster (~ 4—5 Myr), but they would more likely 
have formed BHs. So on one hand the secondary cannot be 
older than '--^ 11 Myr (i?2 = 8 Rq and i?2/-R2,iob = 1) and 
most likely it is not younger than 6 Myr (i?2 ~ 6 Rq 
and i?2/-R2.iob = 0.75). We conclude that the secondary is 
in the late stage of its evolution on the main sequence and 
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probably close to filling its Roche lobe (see also Levine et 
al. 2000). 

We perform the set of calculations for different radii 
of the massive component of LMC X-4 {R2/ R2,\oh = 
0.75 — 0.9) for various efficiencies of tidal interactions 
(-Ftid,rad = 1, 10). The rcsults are shown in Figure [TH The 
orbital decay rate increases with time as the massive com- 
ponent expands along the MS and approaches its Roche 
lobe. The time to reach contact (at which point calcula- 
tions are stopped) decreases with increasing effectiveness 
of tidal forces. For comparison we show the observed or- 
bital decay for LMC X-4, which falls within the model 
with standard tidal interactions efficiency (i^tid,rad = 1)- 
We conclude that in the case of LMC X-4 there is no need 
for the increased efficiency of tidal interactions, and there- 
fore we adopt i^tid.rad = 1 for massive stars with radiative 
envelopes for our standard model value. 

9. X-RAY MODELING 
9.1. X-ray luminosity calculations 

In our study we consider only accreting binaries with 
NS and BH primaries, which are brighter than some X- 
ray luminosity cut Lx,cut- This cut may correspond to a 
detection limit of a particular set of observations. Typ- 
ical -/jx,cut values for most current Chandra observations 
arc in the range 10'^'' — 10'^^ erg s~^. At these high lumi- 
nosities in the Chandra sensitivity range (~ 0.3 — 7 keV) 
the only WD accretors will be supersoft sources, which 
are easily identifiable from their X-ray spectra and arc 
thought to have most of their X-ray emission coming from 
nuclear burning rather than gravitational energy release 
(see Kuulkers et al. 2003 for a review of the X-ray proper- 
ties of WD accretors). Although, for some deep Galactic 
exposures Chandra has reached levels of ~ lO'^" erg s~^ 
(e.g., Galactic Center image of Muno et al. 2003) and a 
contribution from cataclysmic variables may also become 
important. The calculation of X-ray luminosities of sys- 
tems with WD accretors is described in a separate study 
(Ruiter, Belczynski & Harrison 2006). 

Binary companions to NS /BHs may lose material either 
through a stellar wind or via RLOF. In the latter case, 
the donors transfer all the material toward the accrctor, 
whereas for the wind-fed systems only a fraction of the 
material is captured by the compact object. We calculate 
the bolometric luminosity (Lboi) arising from the accretion 
onto a compact object. The accretion rate is based on the 
secular averaged mass accretion rate. If a system is de- 
tached then we use the wind mass accretion rate (eq. I35p , 
and if system is semi-detached the RLOF accretion rate 
is used (eq. [5(11) . Wc do not calculate X-ray luminosities 
arising from the accretion in dynamically unstable phases, 
since the timescales are very short and additionally X-ray 
emission would be highly absorbed due to large optical 
depths in the CE. The Lboi is calculated from 
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where the radius of the accrctor is 10 km for a NS and 
3 Schwarzschild radii for a BH, and e gives a conversion 
efficiency of gravitational binding energy to radiation asso- 
ciated with accretion onto a NS (surface accretion e = 1.0) 
and onto a BH (disk accretion e = 0.5). 



For RLOF-fcd systems wc make a distinction between 
persistent and transient X-ray sources. All wind-fed sys- 
tems are considered as persistent X-ray sources. The issue 
of the wind-fed XRBs with massive Be companions and 
their outburst behavior is discussed in § 9.2. 

RLOF-fed systems are subject to a thermal disk instabil- 
ity and may appear either as persistent or transient X-ray 
sources depending on the mass transfer rate. A system be- 
comes a transient X-ray source when the RLOF rate falls 
below a certain critical value Mdisk- We use the work of 
Dubus et al. (1999) for H-rich disks (see their cq.30) and 
the study of Menou, Perna, & Hernquist (2002) for disks 
with heavier elements (see their eqs.1-4) 
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where A'/acc is accretor mass in Mq, i?disk is a maximum 
disk radius (2/3 of accretor Roche lobe radius) in 10^" cm. 
Constants are: Ci = C j (5 x 10~^), with C being radiation 
parameter of typical value 5 x 10~^; ao.i = a/0.1, with 
a being a viscosity parameter. Following Menou et al. 
(1999) we adopt a = 0.1 for all types of donors since there 
is empirical evidence from dwarf nova outbursts that this 
is the right order of magnitude for the viscosity parame- 
ter. The same value of a is used to derive the critical mass 
transfer rate for H-rich disks (Dubus et al. 1999). H-rich 
donors are the stars with types A'l = 0,1,2,3,4,5,6,16, 
He-rich donors are K-, = 7, 8, 9, 10, 17, CO-rich donors are 
K\ = 11, while we apply formulas for 0-rich type donors 
to ONe WDs {K; = 12). 

We adopt a semi-empirical approach to calculate quies- 
cent X-ray luminosities of transient NS RLOF-fed sources, 
since little is known about the emission mechanism during 
quiescence. It is not certain if the emission arises from a 
low level accretion or a deep crustal heating (for a detailed 
discussion see Belczynski & Taam 2004b, and references 
therein). Using X-ray studies of Galactic transient systems 
with NS accretors (e.g., Tavani & Arons 1997; Rutledge 
et al. 2001; Campana & Stella 2003; Jonker, Wijnands & 
van der Khs 2004; Tomsick et al. 2004; Campana 2004), 
we adopt lO'^^ erg s~^ as a lower limit for the hard X-ray 
luminosity, above 2 kcV. However, it was shown that the 
average luminosity level can be higher ^ 10'^^ erg s"""^ (e.g., 
Rutledge et al. 2002; Jonker et al. 2004). We adopt an X- 
ray luminosity level of 10^^ — 10'^^ erg s~^ above 2 keV. 
Furthermore, we assume that the quiescent NS transient 
X-ray luminosities are evenly distributed (in logLx) in the 
above range. 

The quiescent emission from BH transient systems is 
likely related to a low level of mass accretion. Recent ob- 
servations of BH transients in their quiescent states (Tom- 
sick et al. 2003) reveal rather hard spectra that are not well 
described by a black body. The observed luminosities are 
found in the range 10"^" — 10'^'^ erg s~^ with a median lu- 
minosity ~ 2 X lO'^^ erg s^^. For BH systems wc also use a 
semi-empirical approach, and wc assume that most (80%) 
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of the quiescent BH transient X-ray luminosities above 2 
keV are evenly distributed in the 10^" — 10"^^ erg s^^ range, 
while the rest (20%) of the systems arc slightly brighter: 
luminosities evenly distributed in the ~ 10"^^ — lO'^'^ erg s^^ 
range (see Fig. 3 of Tomsick et al. 2003). Both of the above 
distributions are uniform in logLx- There are some indi- 
cations that the highest quiescent luminosities are found 
in the longest period systems (e.g. Garcia et al. 2001), but 
we do not implement this effect until confirmed by more 
observations. 

RLOF-fed transient systems in outburst reach high 
(close to Eddington) X-ray luminosities. We introduce a 
factor r]out describing the fraction of the critical Eddington 
luminosity a given system has reached. The long period 
systems, with orbits that are sufficiently extensive for a 
large accretion disk to be formed, are usually found to 
emit at the Eddington luminosity (iedd) during outburst, 
while the outburst luminosities of short period systems are 
lower by about an order of magnitude. The correction fac- 
tor to an X-ray luminosity at outburst corresponding to 
'7out = 0.1 and 77out = 1 for the short and long period sys- 
tems is applied respectively. The critical periods, above 
which the Eddington luminosity is adopted, are taken to 
be 1 day and 10 hrs for NS and BH transients in outburst, 
respectively (Chen, Shrader & Livio 1997; Garcia et al. 
2003; see also appendix Al in Portegies Zwart, Dewi & 
Maccarone 2004). 

In order to decide if a given transient system is in an ac- 
tive (outburst) state or inactive (quiescent) state the disk 
duty cycle (DCdisk; the fraction of a time a given system 
spends in the outburst) must be known. However, the 
disk instability theory cannot provide a reliable estimate 
of DCdisk- Empirically it is thought that DCdisk < 1% 
(e.g., Taam, King & Ritter 2000). We adopt DCdisk = 1% 
(probability of finding a system in outburst) in our cal- 
culations and use Monte Carlo to decide the state of a 
transient system. In practice when we study a stellar pop- 
ulation the information for all X-ray binaries is extracted 
at some specified time (time slice). Once a given system 
is identified as a transient (see eq. I84p a random number 
(flat probability distribution) is drawn from the range 0- 
1. If the number is smaller than 0.01 (1% probability) the 
system is then in outburst, otherwise it is in quiescence. 
The appropriate X-ray luminosity is then assigned to the 
system (see eg. I86|) . Alternatively, we use a phenomeno- 
logical model for the duty cycle developed by Portegies 
Zwart et al. 2004. The model is based on the observations 
available for the Galactic BH transient systems. In partic- 
ular comparison of the recurrence time and the decay time 
combined with the observed peak outburst energy allows 
to calculate the time in which system is brighter than a cer- 
tain critical X-ray luminosity. Speciflc application of that 
model will be discussed in the forthcoming paper on the 
evolution of X-ray luminosity function in starburst galax- 
ies (Belczynski et al. 2006, in preparation). 

Finally, the bolometric accretion luminosity is converted 
to an X-ray luminosity in a specific energy range. We 
perform the conversion to the 0.3 - 7 keV range, which 
may be used directly for comparison with Chandra ob- 
servations. For all the persistent RLOF-fed sources, all 
wind-fed sources and the transients in the outburst stage. 
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where accretion is the dominant contributor to the ob- 
served luminosity, we apply a bolometric correction (jyboi)- 
For all quiescent transients the bolometric correction is 
not needed since we adopted their X-ray luminosities di- 
rectly from observations. For different types of systems we 
estimate the correction to be: 
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Corrections were obtained from: La Barbera et al. (2001) 
for wind-fed NS systems; from Di Salvo et al. (2002) and 
Maccarone & Coppi (2003) for RLOF-fed NS systems; and 
from Miller et al. (2001) for BH systems. These bolomet- 
ric corrections will be applicable for the typical Chandra 
observations of external galaxies. For deeper observations, 
where the lower luminosity cutoffs are below a few percent 
of the Eddington limit, the objects make spectral state 
transitions (see Maccarone 2003 and references within), 
and the bolometric corrections are much largei0. 

Combining all of the above information, we can calcu- 
late the X-ray luminosity of synthetic X-ray binaries from 

10^^ — 10^^ all quiescent NS transients 
10^" - 10^2 80% quiescent BH transients 
= < 10^^ - 10^3 20% quiescent BH transients 
f?boifyout-^cdd outburst NS/BH transients 
??boi-^boi persistent (RLOF and wind) 

(86) 

where is expressed in erg s ^ and Ledd represents the 
Eddington luminosity. Note that the X-ray luminosity is 
calculated directly from the mass transfer rate only for per- 
sistent sources. On the other hand, we adopt the above 
empirical description for transient sources since the rela- 
tion between the quiescent, outburst, bolometric luminosi- 
ties, and duty cycle are uncertain due to the mass loss from 
the system during the outburst state. Evidence for such 
mass loss in the form of jets and/or wind have been ob- 
served in, for example, a Galactic BH transient GRS 1915- 
105 (Dhawan, Muno & Remillard 2005; Truss & Done 2006 
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9.2. High Mass X-ray Binaries: Be Star Transients 
9.2.1. Observational Overview 



High mass X-ray binaries consist of a compact object 
(either a NS or a BH) orbiting a massive star. Both galac- 
tic and extra-galactic populations of HMXBs are known 
(Liu, van Paradijs & van den Heuvel 2000, 2005). The 
majority of HMXBs (about 2/3; see Liu et al. 2000, 2005; 
Hayasaki & Okazaki 2005) are so-called Be/X-ray binaries, 
in which the primary is a Be star, orbiting a magnetized 
NS. Orbits are generally wide with a moderate eccentric- 
ity. The compact star accretes from the wind of a massive 
main sequence or subgiant Be (spectral types B3-0 with 
Balmer emission lines; Zorec & Briot 1997) companion. 
Many of these systems show transient behavior (see be- 
low). The remaining HMXBs are those in which the pri- 
mary is a supergiant, so called SG/X-ray binaries (e.g., Liu 

^However, note that the quiescent X-ray luminosities are not affected since they are adopted directly from the deep observations. 
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et al. 2000). For these systems the compact object either 
accretes from the wind of the supergiant, or in brighter 
systems through RLOF (possibly atmospheric but not al- 
ways) via an accretion disk. 

If indeed some HMXBs are confirmed to be evolving 
through stable RLOF, it should pose a useful constraint 
on the development of a delayed dynamical instability. In 
general, it is expected that mass transfer from a much 
more massive donor to a low mass companion is dynam- 
ically unstable and leads to the formation of a CE (see 
§ 5) that ends HMXB phase. It has been shown that if 
a H-rich donor is ~ 3 times more massive than a com- 
pact star accretor (see §5.1) the RLOF will lead to CE 
phase. For adopted the maximum NS mass adopted here 
(2.5 Mq) we predict that only stars of spectral type later 
than B3 (masses smaller than 7.5 M0) could be in dynam- 
ically stable RLOF with NS accretors. If a higher mass 
donor is found in a HMXB with a solid case for ongoing 
RLOF, then either (i) compact object mass is higher (e.g. 
BH) , or (ii) the system is in the phase of short-lived atmo- 
spheric RLOF and will soon end up in CE phase, or (Hi) 
the understanding of development of dynamical instabil- 
ity is incomplete and the observations could be used to set 
new limits. 

Some Be/X-ray binaries (Be XRBs) are persistent 
sources (varying by less than a factor of ~ 10) observed at 
low luminosity levels ~ 10'^^ — 10'^^ erg sec""'^ (e.g.. Van 
Bever & Vanbeveren 2000; Okazaki & Negueruela 2001). 
However, most Be XRBs show periodic outbursts and are 
called transient Be XRBs. Transient Be XRBs exhibit 
two different types of outbursts (e.g., Bildsten et al. 1997; 
Okazaki & Negueruela 2001; Hayasaki & Okazaki 2005; 
Baykal et al. 2005): 

- Type I (normal) outbursts, which are of moderate inten- 
sity (Lx ^ lO'^^ — 10^^ erg sec~^) and they appear to be 
related to the orbital period. It is generally accepted that 
these outbursts are associated with the periastron passage 
of a NS, and are explained by the increased accretion from 
the Be star wind at periastron. 

- Type II (giant) outbursts, with Imninosities reaching 

<; 10"^^ erg sec~^, are irregular, and although they 
seem to appear shortly after the periastron passage, they 
do not exhibit any other correlations with the orbital pe- 
riod. Although the origin of the Type II outbursts remains 
unknown, it was suggested that the outflow from the Be 
star may lead to the formation of a transient accretion disk 
around the NS. Disk accretion results in higher X-ray lu- 
minosities than direct surface wind accretion (see Bildsten 
et al. 1997 for a discussion and references). Some systems 
show both types of outbursts, e.g., A 0535-1-262 (Motch 
et al. 1991; Finger, Wilson & Harmon 1996), V0332-H53 
(Stella, White & Rosner 1986) or 4U 0115+634 (Baykal et 
al. 2005). 

9.2.2. Modeling 

Type I outbursts are averaged out of our calculations 
if we use the orbit-averaged wind accretion model (see 
§4.2). In the general (arbitrary eccentricity) wind accre- 
tion model (see §4.1) Type I outbursts are a natural out- 
come. However, it was noted (Avni & Goldman 1980) that 
the transient phenomenon may be difficult to explain. 

We construct a simple phcnomenological model for Type 
II outbursts in order to be able to assess the influence of 



this transient activity on XRB population characteristics. 
For a system to be a potential Type II Be XRB outburster 
we require: 

- a binary with a NS or a BH accretor and a massive MS 
(Ki = 1) or subgiant (A'l = 2) donor (Af > 8 M©, spectral 
type earlier than B3), 

- that the system is tight enough so it appears as a HMXB 
with a persistent (outside outbursts) wind accretion lead- 
ing to an X-ray luminosity greater than Lx,bg- We allow 
Lx,Be to change within the range 10^^ — 10'^^ erg sec~^. 

Furthermore, only a fraction (/se) of donors in the above 
binaries are Be stars (as opposed to a regular B stars), and 
can potentially trigger the Type II outbursts. To provide 
an upper limit on the contribution of bursting HMXBs to 
the XRB population one may choose /bo = 1- For detail 
studies, the value of /bo may be constrained based on the 
age of a massive star (McSwain & Gies 2005) or its spec- 
tral type and luminosity class (Zorec & Briot 1997). Since 
little is known about the duty cycle of Type II outbursts, 
we allow the duty cycle to change within a wide range 
DCbc — 0.1 — 0.5 and use Monte Carlo to decide whether 
the system is in outburst or in quiescence. Here, DCbb 
gives the fraction of a time a given system spends in the 
outburst. An orbit averaged X-ray luminosity (direct wind 
accretion) is used for quiescence (Tyboi = 0.15, 0.8 § 9.1), al- 
though thermal emission from a NS is also observed in 
some systems. For systems in the Type II outburst the 
X-ray luminosity is taken to be uniformly distributed in 
the range = 10'^'' — lO'^^ erg sec~^. We adopt bolo- 
metric correction factors: ?7boi = 0.15,0.8 for NS and BH 
accretors, respectively (see §9.1). 

The X-ray modeling will be further developed as we pro- 
ceed with the studies of the Galactic and extragalactic 
X-ray binary populations (e.g., Bclczynski et al. 2006, in 
preparation) . 

10. SUMMARY 

We have presented a detailed description of the updated 
StarTrack population synthesis code. The code is be- 
ing used to study populations of different varieties of bi- 
naries hosting compact objects. The code has been cal- 
ibrated and tested against different sets of observations 
and detailed evolutionary calculations and the results are 
presented here. The updated version of StarTrack was 
already used in several studies of compact object bina- 
ries and XRBs. StarTrack allows for evolution of stel- 
lar systems with a wide variety of different initial condi- 
tions (IMF, metallicity, star formation history) and for a 
number of different evolutionary models, subject to the 
parametrization of the input physics. 

The StarTrack code can be compared to the BSE pop- 
ulation synthesis code (Hurley et al. 2002). StarTrack 
incorporates the same single star evolutionary formulas 
(Hurley et al. 2000) as the BSE code, however we extend 
the original formulas to (i) include wind mass loss rates 
from low- and intermediate-mass main sequence stars (for- 
mation of prc-LMXBs, Bclczynski & Taam 2004b); (ii) ac- 
count for the late evolution of low-mass helium stars (im- 
portant for formation of double neutron star systems, see 
Bclczynski et al. 2007 and references therein); and (Hi) 
calculate final masses of neutron stars and black holes, 
based on recent hydrodynamical calculations (e.g., Bel- 
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czynski et al. 2004b). For the treatment of tidal inter- 
actions we use the same equations (ODEs) as in Hurley's 
code, but we employ a fifth-order Runge-Kutta scheme 
with truncation error monitoring and adaptive step-size 
control integration of the ODEs instead of simple multi- 
plication of the derivatives by the evolutionary timestep 
(Euler method). We also adopt convective tides that are 
more efficient (by a factor of 10-100) based on the ob- 
servational calibration discussed in section 8.2. This will 
have a significant effect on the evolution of close binaries 
with low mass (convective) stars. The calculation of X-ray 
luminosities for transient systems with NS and BH accre- 
tors is much more comprehensive in StarTrack. We use 
both the recent theoretical work and observations of low- 
and high-luminosity X-ray sources, to calibrate and test 
our approach (Belczynski & Taam 2004b; Belczynski et 
al. 2004a). The compact object masses formed in core 
collapse are calculated differently, and in particular we ac- 
count for possibility of direct BH formation (no natal kick, 
no mass loss), with maximum BH masses formed reaching 
~ 10 — 20 Mq depending on metallicity and adopted wind 
mass loss rates (e.g., Belczynski et al. 2002; Belczynski et 
al. 2004b), a result that is consistent with maximum BH 
mass estimates in Galactic BH binaries (e.g., ~ 15 Mq for 
GRS 1915; - 10 - 19 M© for Cyg X-1; see Orosz 2003) In 
contrast, in the BSE code all compact objects (including 
BH) have masses below ^2.5 Mq for the entire spectrum 
of initial progenitor masses and different metallicities (see 
Fig. 20 of Hurley et al. 2000), a result that cannot be recon- 
ciled with the current estimates of BH masses. Also more 
recent (Arzoumanian et al. 2002; Hobbs et al. 2005) natal 
kick distributions are used here as compared to Lyne & 
Lorimer (1994) in BSE. The above will affect the post-SNa 
binary orbit, and subsequent evolution of massive binaries. 
For example, the effect of natal kicks on population of dou- 
ble compact objects is rather dramatic and was quantified 
in Belczynski et al. (2002c). The treatment of nuclear mass 
transfer rate is different in the two codes. In BSE it is cal- 
culated using a formula that keeps the donor star within 
its Roche lobe. The formula is calibrated to keep the mass 
transfer steady. In StarTrack, we use the radius- mass 
exponents for the donor and its Roche lobe along with 
an estimate of the evolutionary donor radius change with 
time to calculate the mass transfer rate (see Sec. 5.1). The 
calibration of the BSE prescription is not discussed in de- 
tail by Hurley et al. (2002). Ruiter et al. (2006) find that 
for intermediate polars (low mass main sequence donors 
with WD accretors) the BSE code results in mass trans- 
fer rates of about 2 orders magnitude lower (as calculated 
by Liu & Li 2006) than the rates predicted by StarTrack 
and the observations of intermediate polars may indicate 
(Muno et al. 2006). For low mass binaries, we use a dif- 
ferent (less efficient) magnetic braking law in our standard 
model. As a result, binary orbits in our model will tend 
to take longer time to decay and initiate mass transfer, as 
compared to BSE models. The StarTrack prescription of 



mass accumulation onto white dwarfs is quite unique (see 
§5.7). Related results of calculations for accretion induced 
collapse and NS formation were presented by Belczynski 
& Taam (2004a) and progenitor models of SN la by Bel- 
czynski et al. (2005b). On the other hand, the BSE code is 
more fitted to work with dynamical codes, following in de- 
tail merger products (and their evolution) of various types 
of binary components. Also, the BSE code is much faster 
than the StarTrack code, and therefore may be used for 
simulations of larger stellar populations. 

In a series of papers that will follow we will address the 
issues of modeling of XRBs, and will focus on the compar- 
ison of synthetic XRB populations with the observed X- 
ray point source populations in nearby galaxies. The code 
is also being used to study populations of binaries with 
NSs and BHs as potential source candidates for ground 
based interferometric gravitational radiation observatories 
(e.g., GEO, LIGO, VIRGO) as well as populations of less- 
massive WD binaries for space-based projects (e.g., LISA). 

Although a number of physical processes governing sin- 
gle and binary evolution remain highly uncertain, the ad- 
vances in observational techniques and new results of mas- 
sive surveys allow now various aspects of stellar evolution 
to be explored. We have incorporated several different evo- 
lutionary models within StarTrack (e.g., different mag- 
netic braking laws or CE prescriptions) making possible 
tests of their validity. For example, one such test may be 
based on a comparison of synthetic and observed X-ray 
luminosity functions for nearby starburst galaxies. 

The StarTrack code described in this paper may be 
used only for the evolution of isolated stars and bina- 
ries, i.e., in stellar systems in which dynamical interactions 
are not important (e.g., field populations, open clusters). 
However, a number of interesting studies may be carried 
out for dense stellar environments, in which both stellar 
evolution and dynamical interactions play an important 
role in the formation of compact object binaries. In par- 
ticular, StarTrack was integrated with a dynamical code 
for these types of studies (for details see Ivanova et al. 
2005). 
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Fig. 1. — Final compact object masses in function of initial mass for single star evolution (with solar metallicity and 
standard winds). Top panel shows the full mass range, and indicates pre-collapse mass of the progenitor star. Bottom 
panel shows the mass range important for NS formation, with different types of remnants marked on the plot. For 
discussion see S 2.3.1. 
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Fig. 2. — Initial rotational velocities of stars used in Star Track calculations. In the top panel we present the fit to the 
observational data from StaufFer & Hartmann (1986). In the bottom panel we show the ratio of the data and the model. 
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Fig. 3. — The diagnostic diagram (top panel) used to decide whether a binary should be evolved on a thermal timescale 
or rather RLOF is dynamically unstable (leading to CE evolution and a potential merger). If the mass ratio at the onset 
of RLOF ((?int) is much greater than the mass ratio at the moment when the orbit starts expanding {q\ow) then the system 
is dynamically unstable, otherwise RLOF on a thermal timescale is assumed. The arrow represents the partial derivative 
of donor radius (equal to the Roche lobe radius) with respect to its mass, and points to the place where the donor is 
expected to regain thermal equilibrium. The bottom panel shows a specific system: a 16 Mq Hertzsprung gap donor with 
a 15 M0 MS companion in an 8-day orbit, for which the diagnostic diagram is plotted. The mass transfer begins on a 
thermal timescale (flat part) and then evolves on a slower nuclear timescale (decline). For more details see §5.2. 
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Fig. 4. — The case of a binary disrupted in a supernova explosion: we present the orbit in the coordinate system /// (for 
details see §6). The line OA is parallel to the vector rif//, while the line OB. the asymptote of the hyperbola, is parallel 
to the vector fif^^. The point O is the focus of the hyperbola. 
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Fig. 5. — RLOF sequence for 16 Mq HG + 15 Mq MS binary. Top panel shows mass transfer rate, middle panel orbital 
period, while bottom panel component mass evolution during the RLOF phase. 
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Fig. 6.^ RLOF sequence for 12 Mq MS + 7.5 M© MS binary. Panels same as in Fig. [3 
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Fig. 7. — RLOF sequence for 10 M0 BH + 5 M0 MS binary. The critical Eddington mass accretion rate onto the BH is 
about 3.1 - 4 X 10"'^ Mq yr~^ Panels same as in Fig. [S] 
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Fig. 8. — RLOF sequence for 7 M0 BH + 2 M0 RG binary. The critical Eddington mass accretion rate onto the BH is 
about 2.2 - 2.6 x IQ-"^ yr^i. Panels same as in Fig. [S] 
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Fig. 9. — RLOF sequence for 1.3 Mq NS + 1.6 M© RG binary. The critical Eddington mass accretion rate onto the NS 
is ~ 1.7 X 10^^ M0 yr^^. Panels same as in Fig. [S] 
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Fig. 10. — RLOF sequence for 1.3 Mq NS + 1 M0 RG binary. The critical Eddington mass accretion rate onto the NS 
is ~ 1.7 X 10^^ M0 yr"^. Panels same as in Fig. [3 
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Fig. 11. — RLOF sequence for 1.4 Mq NS + 2.8 Mq evolved He-star binary. The critical Eddington mass accretion rate 
onto the NS is ~ 2.9 x 10"^ Mq yr"^. Panels same as in Fig. [5] Note the very short duration of this RLOF phase; the 
(finite) timesteps taken by the code may be seen through lines showing orbital period and donor mass. 




12. — RLOF sequence for 1.4 Mq NS + 3.6 Mq evolved He-star binary. The critical Eddington mass accretion rate 
the NS is 2.9 x 10"** M© yr"^ Panels same as in Fig. [S] 
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Fig. 13. — Tidal calibration calculation for the open cluster M67. The figure shows the period-eccentricity plane with 
the population of main sequence binary stars at 3.98 Gyr, the current age of the cluster. Bottom and middle panels 
show the results of evolution with increased tidal interactions (-Ftid.con = 100, 10, respectively) as opposed to the standard 
prescription, presented on the top panel (^tid,con = !)• Note the increase of cutoff period (the longest period circular 
binary in a given sample) with increasing i^tid,con- The observed cutoff period for M67 is Pcut — 10 — 12 days. For more 
details see S 8.2.1. 
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Fig. 14. — Tidal calibration calculation for the high-mass X-ray binary LMC X-4. The observed orbital decay rate for 
LMC X-4 is —9.8 x IQ-^yr-^ (marked with dotted line). Predicted decay rates for different radii of the main sequence 
secondary in respect to its Roche lobe {RijRi^xob = 0.75, 0.8, 0.9) are shown for Ftid.rad = 1, 10- Each track is displayed 
after the massive star reaches the relative radius corresponding to a given track and for all tracks the starting point is 
then set to time=0. For more details see 5 8.2.2. 



